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1  Introduction 


Purpose 

The  purpose  of  this  report  is  to  provide  a  complete  set  of  documentation  for 
the  FEMWATER  groundwater  model.  The  intended  users  of  the  manual  will 
have  a  wide  range  of  technical  experience  and  have  widely  different  needs  that 
the  model  and  documentation  will  fulfill.  While  it  is  impossible  to  address  all 
needs  in  the  body  of  one  document,  this  report  has  been  written  to  provide  useful 
information  for  all  users.  Sophisticated  users  will  find  descriptions  of  the 
numerical  techniques  and  a  complete  set  of  governing  equations  that  form  the 
theoretical  basis  of  the  model.  The  casual  user  will  find  examples  of  several 
types  of  problems  that  it  is  hoped  will  include  their  type  of  application.  With 
little  effort  the  casual  modeler  will  be  able  to  follow  the  examples  provided  and 
spend  little  time  with  problem  setup. 


FEMWATER  Origins 

In  the  early  1990’ s,  the  Athens  Laboratory  of  the  U.S.  Environmental 
Protection  Agency  (AERL)  and  the  U.S.  Army  Engineer  Waterways  Experiment 
Station  (WES)  conducted  independent  evaluations  of  a  wide  variety  of 
groundwater  models  to  determine  which  existing  groundwater  models  could  be 
adopted  for  use  in  their  in-house  applications.  AERL  was  interested  in  adopting 
a  three-dimensional  (3-D)  variably  saturated  model  for  wellhead  protection  use 
that  could  model  irregular  geometries.  WES  was  interested  in  the  same 
capabilities  but  for  the  purposes  of  conducting  groundwater  remediation  studies 
at  contaminated  Department  of  Defense  sites  and  for  salinity  intrusion 
applications  in  U.S.  Army  Corps  of  En^neers  navigation  projects.  The 
independent  evaluations  by  both  agencies  resulted  in  the  selection  of  the 
3DFEMWATER  (Yeh  1987b)  and  3DLEWASTE  (Yeh  1990)  models  for  future 
development  and  implementation  within  their  agencies.  Once  it  became  known 
to  the  agencies  that  they  had  similar  interests  and  research  responsibilities,  a 
cooperative  research  agreement  was  reached  and  work  began  on  a  single 
groundwater  modeling  system  to  support  both  agencies.  FEMWATER  is  the 
name  of  the  developed  model. 


Chapter  1  Introduction 


1 


FEMWATER  is  a  modem  implementation  of  the  two  older  models, 
3DFEMWATER  (flow)  and  3DLEWASTE  (transport).  FEMWATER  was 
formed  by  combining  the  two  codes  into  a  single  coupled  flow  and  transport 
model.  The  3DFEMWATER  and  3DLEWASTE  models  were  originally  written 
by  Dr.  Gour-Tsyh  (George)  Yeh  at  Pennsylvania  State  University  while 
FEMWATER  was  written  as  a  collaborative  effort  between  Dr.  Yeh  and  Dr. 
Hsin-Chi  (Jerry)  Lin  at  WES. 

The  improvements  implemented  in  FEMWATER  are  numerous.  First,  the 
entire  program  stmcture  was  changed  to  allow  its  integration  into  the 
Department  of  Defense  Groundwater  Modeling  System  (GMS).  The  GMS 
contains  a  state-of-the-art  graphical  user  environment  that  allows  efficient  model 
setup  and  visualization  (Engineering  Computer  Graphics  Laboratory  (ECGL) 
1996).  This  was  a  particularly  onerous  task  in  the  older  implementation  of  the 
model  since  it  suffered  from  the  common  limitations  of  older  FORTRAN  codes. 
Second,  a  series  of  new  solvers  were  added  to  replace  the  previously  used  block 
iterative  solver.  The  new  solvers  allow  an  arbitrary  node  numbering  scheme  that 
allows  easier  graphical  user  interface  connections  and  still  enjoy  improved 
computational  efficiency.  Third,  density-driven  (salinity)  transport  capability 
was  added  to  allow  salinity  intrusion  studies  in  coastal  aquifers.  This  required 
the  coupling  of  flow  and  transport  within  a  common  model  which  is  the  last 
major  improvement.  Previous  versions  separated  the  flow  and  transport 
calculations. 


Formulation  of  FEMWATER 


FEMWATER  is  designed  to  solve  the  following  system  of  governing 
equations  along  with  initial  and  boundary  conditions,  which  describe  flow  and 
transport  through  saturated-unsaturated  porous  media.  The  governing  equations 
for  flow  are  basically  the  modified  Richards  equation,  which  is  derived  in 
Appendix  A.  The  equation  is  as  follows: 


Governing  equations  for  flow 


p  dh  ^ 
Po  dt 


K 


P  ^ 
Vh  +—Vz 

Po  ) 


p* 

Pc 


(1) 


F 


-1-  P'd  + 


dh 


where 

F  =  storage  coefficient 


(2) 
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h  =  pressure  head 


t  =  time 


K  =  hydraulic  conductivity  tensor 
z  =  potential  head 
q  =  source  and/or  sink 

p  =  water  density  at  chemical  concentration  C 
Po  =  referenced  water  density  at  zero  chemical  concentration 
p  *  =  density  of  either  the  injection  fluid  or  the  withdrawn  water 
6  =  moisture  content 

a'  =  modified  compressibility  of  the  medium 
P'  =  modified  compressibility  of  the  water 
lie  =  effective  porosity  of  the  medium 
S  =  degree  of  saturation 
The  hydraulic  conductivity  K  is  given  by 


K  = 


(p/po)  Peg.  . 


p/po 


^soK 


(3) 


where 

p  =  dynamic  viscosity  of  water  at  chemical  concentration  C 
Po  =  referenced  dynamic  viscosity  at  zero  chemical  concentration 
k  =  permeability  tensor 
ks  =  saturated  permeability  tensor 
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kr  =  relative  permeability  or  relative  hydraulic  conductivity 

Kso  =  referenced  saturated  hydraulic  conductivity  tensor 

The  referenced  value  is  usually  taken  at  zero  chemical  concentration.  The 
density  and  dynamic  viscosity  of  water  are  functions  of  chemical  concentration 
and  are  assumed  to  take  the  following  form 


Po 


=  aj  +  ajC  +  +  a^C^ 


(4) 


and 


^  =  a5+a6C  +  a7C^+a8C 

H'o 


3 


(5) 


where  a^,  a^, ag  are  the  parameters  used  to  define  concentration  dependence  of 
water  density  and  viscosity  and  C  is  the  chemical  concentration. 

The  Darcy  velocity  is  calculated  as  follows 


V  = 


-K 


— Vh  + Vz 

VPo 


(6) 


Initial  conditions  for  flow  equation.  The  initial  conditions  for  the  flow 
equation  are  given  by  Equation  7: 


h  =  hi(x,y,z)  inR,  (7) 

where  R  is  the  region  of  interest  and  hj  is  the  prescribed  initial  condition,  which 
can  be  obtained  by  either  field  measurements  or  by  solving  the  steady  state 
version  of  Equation  (1). 

Boundary  conditions  for  flow  equation.  The  boundary  conditions  for  the 
flow  equation  are  given  in  the  following  equations. 


a.  Dirichlet  conditions: 


h  =  hd(Xb,yb,Zb,t)  onB^, 


(8) 
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b.  Neumann  conditions: 


•n-K-yVh  =  q„(Xb,yb,Zb,t)  onB„, 


(9) 


c.  Cauchy  conditions: 


-n  K 


— Vh  +  Vz 

'.P  ) 


=  qc(Xb’yb>Zb’0  onB,, 


(10) 


d.  Variable  conditions  during  precipitation  period: 


h  =  h  (Xb,yb,Zb,t)  onB^, 


(11) 


or 


-n  K 


^Vh  +  Vz 


=  qp(Xb.yb»Zb>t)  onBy. 


J 


(12) 


e.  Variable  conditions  during  nonprecipitation  period: 


h  =  h  (Xb,yb,Zb,t)  onB^, 


(13) 


or 


h  =  h„(Xb,yb,Zb,t)  onB^, 


(14) 


or 


-n  K 


^Vh  +  Vz 


=  qe(Xb’yb’Zb,t)  onB,, 


J 


(15) 


where 

(Xb,yb,Zb)  =  spatial  coordinate  on  the  boundary 
n  =  outward  unit  vector  normal  to  the  boundary 
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hj  =  Dirichlet  functional  value 
=  Neumann  flux  value 
=  Cauchy  flux  value 
Bjj  =  Dirichlet  boundaiy 
=  Neumann  boundaiy 
Bj,  =  Cauchy  boundary 
Bv  =  variable  boundary 
hp  =  ponding  depth 

qp  =  throughfall  of  precipitation  on  the  variable  boundary 
hm  =  minimum  pressure  on  the  variable  boundary 
qe  =  evaporation  rate  on  the  variable  boundary 

Only  one  of  Equations  1 1-15  is  used  at  any  point  on  the  variable  boundary  at 
any  time. 


Governing  equations  for  transport 

The  governing  equations  for  transport  are  derived  based  on  the  continuity  of 
mass  and  flux  laws  as  given  in  Appendix  A.  The  major  processes  are  advection, 
dispersion/diffusion,  adsorption,  decay,  biodegradation,  and  source/sink. 


ds 

e—  +  p,—  +  \-\C-V-{GD-'^C)  = 


a  — -i-A 

dt 


(0C+p,s)-(eK,C+p,K,S)- 


m - qC  + 

p 


dh 

F — -t-— V-Vl 
dt  p 


(p^ 

— 

de^ 

IpJ 

dt  ^ 

(16) 


S  =  K^C  for  linear  isotherm 


(17) 
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for  Langmuir  isotherm 


(18) 


.S^KC 
1  +  KC 

S  =  KC"  for  Freundlich  isotherm 


(19) 


where 

6  =  moisture  concentration 

=  bulk  density  of  the  medium  (M/L^) 

C  =  material  concentration  in  aqueous  phase  (M/L  ) 

S  =  material  concentration  in  adsorbed  phase  (M/M) 
t  =  time 

V  =  discharge 

V  =  del  operator 

D  =  dispersion  coefficient  tensor 

a'  =  compressibility  of  the  medium 

h  =  pressure  head 

X  =  decay  constant 

m  =  q  Cin  =  artificial  mass  rate 

q  =  source  rate  of  water 

Cin  =  material  concentration  in  the  source 

Kw  =  first  order  biodegradation  rate  constant  through  dissolved  phase 
Ks  =  first  order  biodegradation  rate  through  adsorbed  phase 
F  =  storage  coefficient 
Kd  =  distribution  coefficient 

Snux  =  maximum  concentration  of  medium  in  the  Langmuir  nonlinear 
isotherm 

n  =  power  index  in  the  Freundlich  nonlinear  isotherm 
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K  =  coefficient  in  the  Langmuir  or  Freundlich  nonlinear  isotherm. 
The  dispersion  coefficient  tensor  D  in  Equation  16  is  given  by 


eD  =  a^|V|5 


w 

+  (aL-aT)-^  +  ajTS 


(20) 


where 

IVI  =  magnitude  of  V 

S  =  Kronecker  delta  tensor 

aj.  =  lateral  dispersivity 

<2^  =  longitudinal  dispersivity 

Om  =  molecular  diffusion  coefficient 

T=  tortuosity 

Initial  conditions  for  transport  equation.  The  initial  conditions  for  the 
transport  equation  are  given  by  Equation  16: 

C  =  Ci(x,y,z)  inR  (2i) 

where  R  is  the  region  of  interest  and  Ci  is  the  prescribed  initial  condition,  which 
can  be  obtained  by  either  field  measurements. 

Boundary  conditions  for  transport  equation.  The  boundary  conditions  for 
the  transport  equation  are  given  in  the  following  equations. 


a.  Dirichlet  conditions: 


C  =  Cd(Xb,yb,Zb)  onBj 


(22) 


b.  Variable  conditions; 

n-(vC-0D-VC)  =  n- VC^(xb,yb,Zb,t)  ifn-V<0  (23) 
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n-(-eD-VC)  =  0  ifn-V>0 


(24) 


c.  Cauchy  conditions: 

n  (VC-eD  Vc)  =  q,(xb,yb,Zb,t)  onB,  (25) 


d.  Neumann  conditions: 


n-(-0D  VC)  =  q„(xb,yb,Zb,t)  onB„ 


(26) 


where 

(Xb,yb,Zb)  =  spatial  coordinate  on  the  boundary 

n  =  outward  unit  vector  normal  to  the  boundary 

Cd  =  concentration  on  the  Dirichlet  boundary 

Cv  =  concentration  of  water  through  the  variable  boundary 

Bd  =  Dirichlet  boundary 

By  =  variable  boundary 

=  total  flux  through  the  Cauchy  boundary  Be 

=  total  gradient  flux  through  the  Neumann  boundaries  Bn 

Since  the  hybrid  Lagrangian-Eulerian  approach  is  used  to  simulate  Equation 
16,  it  is  written  in  the  Lagrangian-Eulerian  form  as 


ie+pA)^=v{eDVc)- 


a  — +  A 

dt 


(ec+p,s)- 


{eK,C+p,K,S)  +  m-^qC+ 


„dh  Po  T-r 

^01 

F— -f-— V-V 

— - 

y  dt  p 

<Po  J 

dt] 

(27) 
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(28) 


v.= 


^  +  Pb^d 


for  linear  isotherm  model 


Dt 


dS  dC  ,  X 

=v(eDVc)- 


^  dh 
a’—^X 
dt 


dC  dt 

(ec+p,s)-(eK,c+p,K,s)+ 


P 

m - qC  + 

P 


dh  p 

r?  *  17  V7 

(p) 

de^ 

r  i - V  *  V 

y  dt  p 

(29) 


V.  =■ 


e 


for  Freundlich  and  Langmuir  models 


(30) 


where  Vd  and  Vf  are  the  retarded  and  fluid  pore  velocities,  respectively;  and 
Dvd()/  D  t  and  DyK)/  D  t  denote  the  material  derivative  of  ( )  with  respect  to  time 
using  the  retarded  and  fluid  pore  velocities,  respectively. 

The  flow  equation.  Equation  1,  subject  to  initial  and  boundary  conditions. 
Equations  8-15,  is  solved  with  the  Galerkin  finite  element  method.  The  transport 
equations.  Equations  27  and  28  or  29  and  30,  subject  to  initial  and  boundary 
conditions.  Equations  21-26,  are  solved  with  the  hybrid  Lagrangian-Eulerian 
finite  element  methods.  DeUiiled  implementation  of  the  numerical 
approximation  of  flow  and  transport  problems  is  given  in  Appendix  B. 


10 


Chapter  1  Introduction 


2  Running  FEMWATER 


File  Organization 

FEMWATER  was  designed  to  be  operated  in  batch  mode.  The  input  for 
FEMWATER  is  organized  into  a  set  of  input  files.  The  output  from 
FEMWATER  is  a  combination  of  screen  and  file  output.  A  summary  of  the 
input  and  output  files  is  shown  in  Table  1  and  Table  2. 


Super  File 

When  FEMWATER  is  launched,  the  user  is  prompted  for  the  name  of  a 
single  input  file.  This  file  is  called  the  “super  file”  and  contains  a  list  of  all  of 
the  appropriate  input  and  output  files  used  in  a  particular  simulation.  Grouping 
the  file  names  in  a  super  file  simplifies  file  management  and  eliminates  the  need 
to  type  the  names  of  all  of  the  files  each  time  a  simulation  is  performed.  The 
format  of  the  super  file  is  shown  in  Figure  1. 


FEMSUP 

/*  File  type  identifier  */ 

GEOM  filename 

/*  Geometry  file  7 

BCFT  filename 

r  Model  file  7 

PRTF  filename 

/*  Printed  output  file  7 

ICHD  filename 

/*  Press,  head  init.  cond.  file  7 

ICMC  filename 

/*  Moist,  cont.(nodal)  init.  cond.  file  7 

ICVL  filename 

/*  Velocity  init.  cond.  file  7 

ICON  filename 

r  Concentration  init.  cond.  file  7 

FLVL  filename 

/*  Velocity  flow  file  (for  trans.  only  sim.)7 

FLPH  filename 

/*  Press,  head  file(for  trans.  Only  sim.)  7 

PSOL  filename 

/*  Press,  head  solution  file  7 

MSOL  filename 

/*  Moist,  cont.(nodal)  solution  file  7 

VSOL  filename 

/*  Velocity  solution  file  7 

CSOL  filename 

/*  Concentration  solution  file  7 

Figure  1.  Super  file  format 


Chapter  2  Running  FEMWATER 


11 


Table  1 

Input  Files 

File  Name 

Description 

Super  File 

T ext  file  containing  a  list  of  all  of  the  input  and  output  files  used  in 
a  FEMWATER  simulation. 

Geometry  File 

Text  file  containing  the  data  describing  the  finite  element  mesh,  i.e., 
nodal  coordinates  and  element  topology. 

Model  File 

Text  containing  analysis  parameters  and  options,  material 
properties ,  boundary  conditions,  and  initial  condition  options. 

Initial  Condition  Files 

Text  or  binary  files  containing  concentration,  pressure  head, 
velocity,  moisture  content  initial  conditions. 

Flow  Files 

Text  or  binary  files  containing  a  previously  computed  flow  solution 
(pressure  head  and  velocity)  which  are  used  to  define  a  3-D  flow 
field  for  a  transport  only  simulation. 

Table  2 

Output  Files 

File  Name 

Description 

Printed  Output 

Text  file  containing  a  summary  of  the  output. 

Pressure  Head 

Text  or  binary  file  containing  the  computed  pressure  heads.  Used 
for  post-processing  or  as  initial  conditions  for  a  subsequent 
analysis. 

Moisture  Content 

Text  or  binary  file  containing  the  computed  moisture  content  at 
nodes.  Used  for  post-processing. 

Velocity 

Text  or  binary  file  containing  the  computed  Darcian  velocities.  Used 
for  post-processing. 

Concentration 

Text  or  binary  file  containing  the  computed  concentrations.  Used 
for  post-processing  or  as  initial  conditions  for  a  subsequent 
analysis. 

The  first  record  in  the  file  is  the  file  type  identifier.  Each  of  the  subsequent 
records  represents  an  input  or  an  output  file.  The  first  field  in  each  record  is  a 
four-character  string  identifying  the  type  of  the  file  listed  in  the  record.  The 
second  field  in  each  record  is  the  name  of  the  corresponding  file.  The  files 
should  be  in  the  same  directory  as  the  super  file. 

Not  all  of  the  files  shown  in  Figure  1  are  required  for  every  simulation. 
Some  of  the  initial  condition  files  are  not  required  depending  on  the  iniHal 
condition  options  specified  in  the  model  file.  Also,  the  user  can  also  specify  in 
the  model  file  not  to  output  some  of  the  solution  files. 


Card  Style  Format 


The  records  used  in  the  super  file  format  shown  in  Figure  1  are 
representative  of  the  formatting  style  used  for  all  of  the  FEMWATER  input  files. 


12 


Chapter  2  Running  FEMWATER 


This  format  is  often  referred  to  as  the  “card  style”  format.  With  this  format,  the 
components  of  the  file  are  grouped  into  logicd  groups  called  “cards.”  Typically, 
each  card  is  a  single  line  or  record;  however,  some  cards  extend  to  multiple 
lines.  The  first  component  of  each  card  is  a  short  name  that  serves  as  the  card 
identifier.  The  remaining  fields  on  the  line  contain  the  information  associated 
with  the  card.  In  some  cases,  such  as  lists,  a  card  can  use  multiple  lines.  All  of 
the  cards  are  assumed  to  be  free-format. 

Several  advantages  are  associated  with  the  card  type  approach  to  formatting 
files: 

a.  Card  identifiers  make  the  file  easier  to  read.  Each  input  line  has  a  label, 
which  helps  to  identify  the  data  on  the  line. 

b.  The  card  names  are  useful  as  text  strings  for  searching  in  a  large  file. 

All  input  lines  of  a  particular  type  can  be  located  quickly  in  a  large  input 
file. 

c.  Cards  allow  the  data  to  be  input  in  any  order  in  many  cases;  i.e.,  the 
order  that  the  cards  appear  in  the  file  is  usually  not  important. 

d.  Cards  make  it  easy  to  modify  a  file  format.  New  data  can  be  included 
simply  by  defining  a  new  card  type.  If  the  new  card  is  optional  (which  is 
typically  the  case  for  new  cards)  old  files  are  still  compatible.  If  an  old 
card  type  is  no  longer  used,  the  card  can  simply  be  ignored  without 
causing  input  errors. 


Other  Files 

Each  of  the  files  listed  in  the  super  file  are  described  in  more  detail  in 
subsequent  chapters.  The  geometry  file  is  described  in  Chapter  3,  the  contents  of 
the  model  file  are  described  in  Chapters  3,  5, 6,  and  7,  the  initial  condition  files 
are  described  in  Chapter  7,  and  the  solution  files  are  described  in  Appendix  C. 
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3  Meshes 


introduction 

The  computational  discretization  utilized  by  FEMWATER  is  a  three- 
dimensional  finite  element  mesh.  The  volumetric  domain  to  be  modeled  by 
FEMWATER  must  be  idealized  and  discretized  into  hexahedra,  prisms,  and/or 
tetrahedra.  Elements  are  typically  grouped  into  zones  representing  different 
stratigraphic  units.  Each  element  is  assigned  a  material  DD  representing  the  zone 
to  which  the  elements  belongs.  When  constructing  a  mesh,  care  should  be  taken 
to  ensure  that  elements  do  not  cross  or  straddle  stratigraphic  boundaries. 


Elements  Supported 


The  types  of  elements  supported  by  FEMWATER  are  shown  in  Figure  2. 
Each  of  the  elements  are  linear;  quadratic  elements  are  not  supported.  Although 
all  three  element  types  are  supported,  tetrahedra  do  not  perform  as  well  as  the 
other  types  and  should  be  avoided  if  possible.  The  numbering  sequence  shown 
in  Figure  2  should  be  used  when  describing  the  elements  in  the  geometry  file. 


Geometry  File  Format 

The  coordinates  of  the  mesh  nodes  and  the  element  topology  are  input  to 
FEMWATER  through  the  geometry  file.  The  format  of  the  geometry  file  is 
shown  in  Figure  3. 
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Figure  2.  The  element  types  supported  by  FEMWATER 


3DFEMGEO 

r  File  type  identifier  */ 

T1  text 

/*  Title,  line  1  ’/ 

T2  text 

r  Title,  line  2  V 

TStext 

/*  Title,  line  3  */ 

GN  id  x  y  z 

r  Nodal  coordinates  */ 

GE8  id  n1  n2  n3  n4  n5  n6  n7  n8  matid 

r  Hex.  Element  */ 

GE6  id  n1  n2  n3  n4  n5  n6  matid 

r  Prism  element  */ 

GE4  id  n1  n2  n3  n4  matid 

/*  Tetrahedral  element  7 

END 

r  End  of  input  data  7 

Figure  3.  Geometry  file  format 

The  cards  used  in  the  geometry  file  are  as  follows: 


Card  Type 

3DFEMGEO 

Description 

File  tvoe  identifier.  Must  be  on  first  line  of  file.  No  fields. 

Repuired 

YES 
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Card  Type 

T1,T2,T3 

Description 

Three  title  lines.  Any  text  can  be  entered  on  these  lines.  The  text  is  used  as  a  banner 
to  the  printed  output. 

Required 

NO 

Format 

T1  text 

T2text 

TStext 

Sample 

T1  Salinity  intrusion  problem 

T2  ACR  Environmental  Sen/ices  Inc. 

T3  June  19, 1995 

Card  Type 

GE8 

Description 

Hexahedral  element.  One  entry  for  each  hexahedral  element  in  the  mesh.  The 
ordering  of  the  nodes  should  correspond  to  the  diagram  in  Figure  2. 

Required 

YES 

Format 

GE8  id  n1  n2  n3  n4  n5  n6  n7  n8  matid 

Sample 

GE8  1 0  847  938  943  928  380  942  835  655  1 

Field 

Variable 

Value 

Description 

id 

NEL 

+ 

The  ID  of  the  node. 

N1 

lE(NEL.I) 

+ 

The  index  of  node  number  1 . 

N2 

IE(NEL,2) 

+ 

The  index  of  node  number  2. 

N3 

IE(NEL.3) 

+ 

The  index  of  node  number  3. 

N4 

IE(NEL,4) 

+ 

The  index  of  node  number  4. 

N5 

IE(NEL.5) 

+ 

The  index  of  node  number  5. 

N6 

IE(NEL.6) 

+ 

The  index  of  node  number  6. 

N7 

IE(NEL,7) 

+ 

The  index  of  node  number  7. 

N8 

IE(NEL,8) 

+ 

The  index  of  node  number  8. 

Matid 

IE(NEL,9) 

The  element  material  index. 

Card  Type 

GE6 

Description 

Wedge  or  prism  element.  One  entry  for  each  wedge  element  in  the  mesh.  The 
ordering  of  the  nodes  should  correspond  to  the  diagram  in  Figure  2. 

Required 

YES 

Format 

GE6  id  n1  n2  n3  n4  n5  n6  matid 

Sample 

GE6  1 0  847  938  943  928  380  942  1 

Field 

Variable 

Value 

Description 

id 

NEL 

+ 

The  ID  of  the  node. 

N1 

IE(NEL,1) 

+ 

The  index  of  node  number  1 . 

N2 

IE(NEL.2) 

+ 

The  index  of  node  number  2. 

N3 

IE(NEL,3) 

+ 

The  index  of  node  number  3. 

N4 

IE(NEL.4) 

+ 

The  index  of  node  number  4. 

N5 

IE(NEL,5) 

+ 

The  index  of  node  number  5. 

N6 

IE(NEL.6) 

+ 

The  index  of  node  number  6. 

Matid 

IE(NEL,9) 

+ 

The  element  material  index. 
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Card  Type 

GE4 

Description 

Tetrahedral  element.  One  entry  for  each  tetrahedral  element  in  the  mesh.  The 
orderina  of  the  nodes  should  correspond  to  the  diaaram  In  Flqure  2. 

Reauired 

YES 

Format 

GE4  id  n1  n2  n3  n4  matid 

Sample 

GE4  10  847  938  943  928  1 

Field 

Variable 

Value 

Description 

id 

NEL 

+ 

The  ID  of  the  node. 

N1 

lEfNEUD 

+ 

The  index  of  node  number  1 . 

N2 

IE(NEL,2) 

+ 

The  index  of  node  number  2. 

N3 

IE(NEL,3) 

+ 

The  index  of  node  number  3. 

N4 

IE(NEL,4) 

+ 

The  index  of  node  number  4. 

Matid 

IE(NEL,9) 

+ 

The  element  material  index. 

Card  Type 

GN 

Description 

Mesh  node.  One  entry  for  each  node  in  the  mesh. 

Repuired 

YES 

Format 

GN  id  X  v  z 

Sample 

GN  83  3482.4  4389.3  34.6 

Field 

Variable 

Value 

Description 

id 

NNP 

+ 

The  ID  of  the  node. 

X 

XfNNP.I) 

± 

The  X  coordinate  of  the  node  fU. 

V 

X(NNP,2) 

± 

The  V  coordinate  of  the  node  ILl. 

Z 

XfNNP,3) 

± 

The  z  coordinate  of  the  node  [LI. 

Card  Type 

END 

Description 

End  of  input  data  marker. 

Repuired 

YES 

Mesh  Generation  Guidelines 

As  with  any  finite  element  problem,  special  care  should  be  used  in  the 
constraction  of  the  mesh  to  avoid  numerical  instability.  Numerical  accuracy  and 
stability  are  often  competing  interests.  Procedures  that  assist  in  achieving  model 
stability  often  do  so  to  the  detriment  of  accuracy.  When  the  need  for 
computational  efficiency  is  added  to  accuracy  and  stability  concerns,  mesh 
generation  can  be  a  difficult  task.  Indeed,  many  consider  it  an  art  form  although 
numerical  efforts  are  under  way  to  minimize  the  artistic  requirements. 

Prior  to  mesh  generation,  the  hydrogeologic  conceptual  model  should  be 
well  defined.  The  stratigraphy  and  associated  variations  in  hydraulic 
conductivity  should  be  weU  understood.  Boundary  conditions,  whether  they  are 
surface  recharge,  salinity  concentrations,  or  wells,  should  be  known  in  detail. 
The  locations  of  remediation  structures  or  other  mitigative  treatments  should  be 
mapped.  Once  all  of  this  information  is  known,  it  is  possible  to  start  the  mesh 
generation  process.  Guidelines  are  now  presented  to  assist  in  the  process. 

The  first  step  in  setting  up  a  mesh  is  to  take  the  solid  model  generated  in  the 
subsurface  conceptualization  phase  and  determine  where  the  external  boundary 


Chapters  Meshes 


17 


conditions  will  be  located.  These  will  essentially  be  water  table  information 
around  the  area  of  interest.  The  external  edges  of  the  model  should  be  located 
where  this  information  is  well  known.  The  next  step  is  to  determine  the  number 
of  elements  and  their  best  distribution  to  solve  the  flow  and  transport  problem. 
Simply  put,  fine  mesh  spacing  should  be  located  where  head  or  concentration 
gradients  are  expected  to  be  maximum.  This  will  certainly  be  near  wells  that 
have  caused  significant  cones  of  depression.  However,  care  must  be  given  to 
gradually  vary  the  size  of  the  elements  to  avoid  numerical  errors.  For  example, 
the  50  percent  rule  should  be  followed  whenever  possible:  the  size  of  an  element 
should  not  differ  from  the  size  of  an  adjacent  element  by  more  than  50  percent. 
Additionally,  every  effort  should  be  made  to  avoid  highly  skewed  or  irregularly 
shaped  elements.  Equilateral  triangles  are  ideally  shaped. 

Once  a  horizontal  (or  surface)  mesh  has  been  developed,  interest  should  shift 
to  the  proper  vertical  element  spacing.  The  same  issues  arise  again  in  the 
vertical  problem.  Fine  resolution  should  be  placed  in  vertical  regions  where 
head  or  concentration  gradients  are  greatest  and  most  particularly  in  the 
unsaturated  zone.  For  example,  if  a  highly  conductive  aquifer  is  adjacent  to  a 
highly  impermeable  aquiclude,  fine  mesh  resolution  is  required  in  the  vicinity  of 
the  interface.  In  general,  there  should  be  a  minimum  of  three  layers  of  elements 
vertically  for  each  distinct  stratigraphic  unit  particularly  if  large  variations  of 
hydraulic  conductivity  occur  in  adjacent  layers.  As  an  aside,  if  large  variations 
in  hydraulic  conductivity  are  required,  no  two  adjacent  layers  should  vary  by 
rnore  than  three  orders  of  magnitude.  If  this  mle  is  violated,  the  solutions  will 
likely  be  inaccurate  and  probably  slow  to  converge.  If  there  are  indeed  sharp 
variations  in  conductivity,  then  many  layers  should  be  used  and  the  values  varied 
gently  over  the  short  distance  in  which  they  change. 

When  setting  up  a  mesh  for  a  transport  analysis,  all  of  the  previous  issues  are 
germane  as  well  as  one  other.  First,  if  hexahedral  elements  are  used,  each 
element  should  be  constracted  such  that  all  of  the  element  faces  are  planar.  This 
is  a  particular  requirement  because  the  particle  tracking  algorithm  used  by 
FEMWATER  may  break  down  if  a  particle  crosses  an  element  face  that  is  not 
planar.  Triangular  faces  are  always  planar,  but  quadrilateral  faces  may  or  may 
not  be  planar.  In  general,  it  is  easiest  to  use  triangular  elements  to  avoid 
problems. 

On  the  subject  of  computational  efficiency,  it  is  important  to  note  that  the 
smallest  number  of  elements  does  not  always  provide  the  fastest  simulation.  It  is 
quite  possible  to  constract  a  model  with  insufficient  numbers  of  elements  to 
characterize  the  problem  adequately,  thereby  creating  a  simulation  that  is  slow  to 
converge.  In  short,  fewer  elements  will  be  used  in  the  calculation,  but  greater 
numbers  of  iterations  will  be  required.  In  general,  if  sufficient  numbers  of 
elements  are  used,  fewer  iterations  will  be  required  to  converge  on  a  solution.  It 
is  better  to  use  large  numbers  of  elements  for  few  iterations  to  get  accurate 
answers  than  to  use  few  elements  for  many  iterations  to  get  inaccurate  and 
possibly  divergent  answers. 
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4  Analysis  Options 


Introduction 

One  of  the  primary  FEMWATER  input  files  is  the  model  file,  which  consists 
of  analysis  options,  material  properties,  boundary,  and  initial  conditions.  The 
first  of  these  groups,  the  analysis  options,  are  described  in  this  Chapter.  The 
material  properties  are  described  in  Chapter  5,  the  boundary  conditions  are 
described  in  Chapter  6,  and  the  initial  conditions  are  described  in  Chapter  7. 


File  Format 

The  set  of  cards  in  the  model  file  corresponding  to  analysis  options  is  shown 
in  Figure  4.  The  file  type  identifier  and  the  title  cards  are  similar  to  the 
corresponding  cards  found  in  the  geometry  file  and  described  in  Chapter  3.  Each 
of  the  remaining  cards  is  described  in  more  detail  in  the  following  sections. 


Run  Option  Parameters 

The  run  option  parameters  designated  on  cards  OPl,  OP2,  OPS,  OP4,  and 
OPS  include  options  for  specifying  the  type  of  simulation,  the  solver,  relaxation 
parameters,  and  sorption  options. 


Type  of  simulation  (0P1) 

The  OPl  card  is  used  to  specify  the  type  of  simulation  to  be  performed  by 
FEMWATER.  The  parameter  KMOD  indicates  the  type  of  simulation  to  be 
conducted. 
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3DFEMWBC 
T1  text 
T2  text 
T3text 
OP1  kmod 

OP2  kssf  ksst  ilump  imid  ipntsf  iquar 
OP3  wf  omef  omit 
OP4  ksorp 
OPS  gg 

IP1  niterf  ncycle  npiterf  tolaf  tolbf 

IP2  nitert  npiterttolbt 

IP3  nitfit  omeftt  epss  epst 

PT1  nxw  nyw  nzw  idetq 

TC1  tmax 

TC2  idt  delt  idtxy 

OC1  ibug  ichng  jopt  kprt/nprint 

OC2  nselt  kpro() 

OC3  ifile  kept  kdsk/npost 

OC4  kselt  ksaveO _ 


/*  File  type  identifier  */ 

/*  Title,  line  1  7 
/*  Title,  line  2  7 
/*  Title,  line  3  7. 

/*  Simulation  type  7 
/'  Solver  options  7 
/*  Weight  factor  7 
/*  Sorption  7 

/*  Upper  bound  for  eigen  value  7 
/’  Iteration  opt.s,  flow  7 
r  Iteration  opt.s,  trans.  7 
/*  Iteration  opts,  coupled  7 
r  Particle  tracking  opts  7 
/*  Total  simulation  time  7 
/*  Time-step  definition  7 
/*  Print  options  7 
/*  Print  options  7 
/*  Format  and  interval  7 
/*  Solution  file  opts  7 _ 


Figure  4.  The  analysis  option  cards  in  the  model  file 


Card  Type 

OP1 

Description 

Type  of  simulation. 

Repuired 

YES 

Format 

OP1  kmod 

Sample 

OP1  10 

Field 

Variable 

Value 

Description 

1 

KMOD 

10 

Flow  simulation  only. 

1 

Transport  simulation  only. 

11 

Coupled  flow  and  transport. 

The  following  options  are  available: 

a.  Perform  a  flow  simulation  only  (KMOD=10). 

b.  Perform  a  transport  simulation  only(KMOD=l).  For  this  case,  a  steady- 
state  or  transient  flow  simulation  must  be  performed  prior  to  the 
transport  simulation.  The  results  of  this  simulation  are  then  input  to 
FEMWATER  as  flow  variables  (pressure  head  or  velocity).  The 
pressure  head  is  required  for  a  transport  simulation.  The  steps  involved 
in  setting  up  the  proper  initial  conditions  are  described  in  Chapter  7. 

c.  Perform  a  coupled  flow  and  transport  simulation  (KMOD=l  1).  With  a 
coupled  flow  and  transport  simulation,  the  user  has  the  option  of 
simulating  either  density-dependent  flow  or  density-independent  flow. 
This  option  is  controlled  by  entering  the  appropriate  parameters  defining 
the  relationship  between  concentration  and  density  and  concentration 
and  viscosity.  These  parameters  are  entered  on  the  MP4  card  described 
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in  the  section  “Concentration  Dependence  Coefficients  (MP4)”  in 
Chapter  5. 

(1)  Density-dependent  flow.  In  this  case,  the  flow  and  transport 
simulations  are  performed  simultaneously.  The  concentration  of  the 
solute  changes  the  density  of  the  solution,  thus  changing  the  flow 
solution.  This  option  should  be  chosen  only  for  density  dependent 
flow  problems  such  as  salinity  intrusion  in  coastal  aquifers.  The 
proper  setting  for  this  case  is  shown  in  sample  problem  5  of  Chapter 
8. 

(2)  Density-independent  flow.  In  this  case,  the  flow  and  transport 
simulation  are  performed  sequentially  for  every  time-step.  The 
proper  setting  for  this  case  is  shown  in  sample  problem  3  of  Chapter 
8. 


Solver  options  (OP2) 

The  OP2  card  is  used  to  select  the  type  of  time  mode  (steady-state  or 
transient)  to  be  used  and  to  set  various  solver  options. 


Card  Type 

OP2 

Description 

Solver  options. 

Reouired 

YES 

Format 

OP2  kssf  ksst  iiump  imid  ipntsf  iquar 

Sample 

OP2 1110  11 

Field 

Variable 

Value 

Description 

1 

KSSF 

0 

Steady-state  flow  simulation. 

1 

Transient  flow  simulation. 

2 

KSST 

0 

Steady-state  transport  simuiation. 

1 

T ransient  transport  simuiation. 

(Note:  KSSF  and  KSST  must  be  set  to  the  same  value) 

3 

iLUMP 

0 

No  mass  lumping. 

1 

Mass  iumpinq. 

4 

IMiD 

0 

No  mid-difference. 

1 

Mid-difference. 

5 

IPNTSF 

1 

The  pointwise  iterative  matrix  solver. 

2 

P.C.G.  method  (polynomial). 

3 

P.C.G.  method  (incomplete  Cholesky). 

6 

IQUAR 

11 

Nodal/nodal  quadrature. 

12 

Nodal/  gaussian  quadrature. 

21 

Gaussian/nodal  quadrature. 

22 

Gaussian/qaussian  quadrature. 

Steady-state  versus  transient.  MWATER  can  be  run  in  either  a  steady- 
state  or  transient  mode.  The  steady  state  mode  is  allowed  only  when  the  flow- 
simulation-only  option  has  been  selected  with  the  OPl  card.  FEMWATER  uses 
the  Lagrangian-Eulerian  finite  element  method  to  solve  the  transport  equation. 
Therefore,  the  steady-state  mode  of  transport  simulation  is  not  allowed  in  this 
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option.  The  transient  mode  must  be  used  when  a  transport  simulation  is  being 
performed. 

Mass  lumping  (BLUMP).  This  parameter  indicates  whether  or  not  the  mass 
matrix  is  to  be  lumped.  With  lumping  (ILUMP=1),  the  solution  is  less  accurate 
but  potentially  more  stable.  For  saturated-unsaturated  flow  computations,  or  if 
negative  concentrations  or  oscillating  solutions  occur,  this  parameter  should  be 
set  to  1.  If  the  computations  are  quite  stable,  particularly  in  largely  saturated 
flow  simulations,  the  parameter  should  be  set  to  zero. 

Mid-difference  (IMID).  This  parameter  indicates  if  the  mid-difference 
method  should  be  used  in  both  the  flow  and  transport  computations.  If  IMID  =1, 
the  mid-difference  method  is  used.  Setting  IMID=1  is  reserved  for  research 
purposes  so  IMID=0  is  the  preferred  setting. 

Solver  Selection  (IPNTSF).  The  following  three  solvers  are  provided  in 
FEMWATER: 

a.  Pointwise  iterative  matrix  solver.  The  pointwise  iterative  matrix  solver 
employs  the  basic  successive  iterative  method  to  solve  the  matrix 
equation,  including  the  Gauss-Seidel  method,  successive 
underrelaxation,  and  successive  overrelaxation.  When  the  resulting 
matrix  is  diagonally  dominant,  the  pointwise  iterative  solver  provides  a 
convergent  solution.  This  solver  is  preferred  because  it  is  more  robust 
than  the  other  two  solvers.  However,  when  the  speed  of  convergence  is 
too  slow,  one  may  wish  to  choose  one  of  the  other  two  solvers. 

b.  Preconditioned  conjugate  gradient  method  (polynomial).  This  solver 
employs  the  conjugate  gradient  method  to  solve  the  matrix  equation.  It 
uses  a  polynomial  as  a  preconditioner.  This  matrix  solver  provides  a 
convergent  solution  when  the  resulting  matrix  is  symmetric  positive 
definite  (SPD).  Theoretically,  the  convergence  speed  is  faster  than  the 
pointwise  iterative  solver.  This  solver  should  be  used  only  when  the 
pointwise  iterative  solver  is  too  slow. 

c.  Preconditioned  conjugate  gradient  method  (incomplete  Choleski).  This 
solver  employs  the  conjugate  gradient  method  using  the  incomplete 
Choleski  decomposition  as  the  preconditioner.  A  convergent  solution  is 
provided  when  the  matrix  is  SPD.  However,  when  the  matrix  is  slightly 
non-symmetric,  the  solver  could  also  give  convergent  solutions.  Its 
speed  of  convergence  is  theoretically  faster  than  the  pointwise  iterative 
solver  and  is  comparable  to  the  polynomial  preconditioned  conjugate 
gradient  method.  This  solver  should  be  used  only  when  the  pointwise 
iterative  solver  is  too  slow.  This  solver  is  generally  but  not  always 
preferred  over  the  polynomial  preconditioned  conjugate  gradient 
method. 
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Quadrature  selection  (IQUAR).  This  parameter  is  an  indicator  of  the  type 
of  quadrature  used  in  the  numerical  integration.  The  following  four  quadrature 
options  are  provided: 

a.  Nodal/nodal  quadramre.  Nodal  quadrature  is  used  for  surface  and 
element  integration. 

b.  Nodal/gaussian  quadrature.  Nodal  quadrature  is  used  for  surface 
integration,  and  gaussian  quadrature  is  used  for  element  integration. 

c.  Gaussian/nodal  quadrature.  Gaussian  quadrature  is  used  for  surface 
integration,  and  nodal  quadrature  is  used  for  element  integration. 

d.  Gaussian/gaussian  quadrature.  Gaussian  quadrature  is  used  for  both 
surface  and  element  integration. 

Gaussian/gaussian  quadrature  yields  the  most  accurate  solution  and  should 
be  used  as  the  default  v^ue.  However,  this  option  may  provide  oscillations  or 
divergence  in  highly  nonlinear  problems.  When  this  occurs,  the  user  should  try 
to  use  the  nodal/nodal  quadrature.  Once  these  options  have  been  used 
unsuccessfully,  the  remaining  options  can  be  tried. 


Weighting  factor  options  (OPS) 

The  OPS  card  is  used  to  select  the  type  of  time  derivative  and  relaxation 
weighting  factors  to  be  used  and  to  set  several  parameters  associated  with  the 
weighting  factor. 


Card  Type 

OP3 

Description 

Weighting  factor  options. 

Repuired 

YES 

Format 

OPS  wf  omef  omif 

Sample 

OPS  1.0  1.0  1.0 

Field 

Variable 

Value 

Description 

1 

WF 

0.5 

Crank-Nicolson. 

1.0 

Backward  difference. 

2 

OMEF 

Iteration  param.  For  nonlinear  flow  and  transport. 

0.0-1 .0 

Underrelaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Overrelaxation. 

3 

OMIF 

Iteration  param.  For  linearized  flow  and  transport. 

0.0-1 .0 

Underreiaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Overrelaxation. 

Weighting  factor  type  (WF).  This  parameter  determines  how  one  would 
evaluate  the  time  derivative  terms  associated  with  the  velocity  in  the  flow 
equation.  Two  types  of  weighting  factors  are  available: 
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a.  Crank-Nicolson  central  (WF=0.5).  When  new  time  derivatives  are 
determined  by  averaging  the  previous  time  derivative  and  an  estimated 
time  derivative,  the  process  is  called  Crank-Nicolson  central  weighting. 

b.  Backward  difference  (WF=  1.0).  When  the  time  derivatives  are 
evaluated  only  at  the  new  time,  the  process  is  called  backward  difference 
weighting. 

A  value  of  WP  equal  to  1.0  (an  implicit  numerical  scheme)  should  be  used 
for  most  practical  problems.  Setting  WT  equal  to  0.5  is  normally  done  for 
research  purposes  to  assess  the  accuracy  of  the  Crank-Nicolson  scheme. 

Relaxation  parameter  for  solving  nonlinear  flow  and  transport 
equations  (OMEF).  WTien  the  flow  and  transport  equations  are  nonlinear,  an 
estimate  of  the  pressure  head  and  the  concentration  is  needed  to  compose  the 
matrix  equation.  There  are  three  options  to  estimate  the  pressure  head  and 
concentration  based  on  previous  guesses  and  newly  obtained  values: 
underrelaxation,  exact  relaxation,  and  overrelaxation.  OMEF  is  a  weighting 
factor  that  is  applied  to  the  newly  obtained  values,  and  a  weighting  factor  of  1.0 
minus  OMEF  is  applied  to  the  previous  guesses.  For  underrelaxation,  a  value  of 
OMEF  between  0.0  and  1 .0  is  used  for  the  newly  obtained  values.  For  exact 
relaxation  OMEF  is  set  equal  to  1.0  and  the  newly  obtained  values  are  used  as 
the  new  guesses.  For  overrelaxation,  a  value  of  OMEF  between  1.0  and  2.0  is 
used  for  the  newly  obtained  values. 

Normally  OMEF  should  be  set  to  1.0.  If  the  convergence  history  shows 
signs  of  oscillation,  then  a  value  of  0.5  should  be  used.  If  the  convergence 
history  shows  a  monotonic  decrease  but  at  a  very  slow  rate,  then  it  should  be  set 
to  between  1.7  and  1.9. 

Relaxation  parameter  for  solving  linearized  flow  and  transport 
equations  (OMIF).  In  order  to  solve  the  linearized  matrix  equations  using  the 
iteration  method,  an  estimate  of  the  solution  is  needed  prior  to  taking  the  next 
iteration.  There  are  three  options  to  estimate  the  solution  based  on  previous 
guesses  and  the  newly  obtained  solution:  underrelaxation,  exact  relaxation,  and 
overrelaxation.  This  is  accomplished  with  an  OMIF  weighting  factor  that  is 
similar  to  the  OMEF  factor  described  in  the  previous  section. 

Normally  OMIF  should  be  set  to  1.0.  If  the  convergence  history  shows  signs 
of  oscillation,  then  a  value  of  0.5  should  be  used.  If  the  convergence  history 
shows  a  monotonic  decrease  but  at  a  very  slow  rate,  then  it  should  be  set  to 
between  1.7  and  1.9. 


Sorption  options  (OP4) 

The  OP4  card  is  used  to  designate  which  model  wiU  be  used  for  the  sorption 
isotherm.  The  selection  of  the  sorption  model  should  be  dictated  by 
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experimental  evidence,  and  it  depends  highly  on  the  type  of  chemicals  and 
subsurface  media. 


Card  Type 

OP4 

Description 

SoroBon  options. 

Repuired 

YES 

Format 

OP4  ksorp 

Sample 

OP41 

Field 

Variable 

Value 

Description 

1 

KSORP 

1 

Linear  isotherm. 

2 

Freundlich  isotherm. 

3 

Lanqmuir  isotherm. 

The  following  three  models  are  available  for  modeling  the  sorption  isotherm: 

a.  Linear.  A  linear  isotherm  is  used  for  the  adsorption  model.  For  salinity 
intrusion  simulations,  a  linear  model  is  sufficient. 

b.  Freundlich.  A  nonlinear  isotherm  (Freundlich  isotherm)  is  used  for  the 
adsorption  model. 

c.  Langmuir.  A  nonlinear  isotherm  (Langmuir  isotherm)  is  used  for  the 
adsorption  model. 

Although  the  Freundlich  isotherm  option  can  be  used  to  simulate  a  linear 
isotherm  by  setting  the  value  of  the  exponent  (n  =  1),  it  is  recommended  that  the 
linear  isotherm  be  simulated  by  using  only  the  linear  isotherm  option.  This  is 
because  the  linear  isotherm  option  makes  use  of  retarded  seepage  velocities, 
which  result  in  a  more  accurate  solution  for  the  particle  tracking  scheme  than  the 
pore  velocities  used  in  conjunction  with  the  nonlinear  adsorption  models. 


Preconditioned  conjugate  gradient  method  (OPS) 

The  OP5  card  is  used  to  provide  an  estimator  for  GG  in  the  preconditioned 
conjugate  gradient  method.  If  the  OPS  card  is  not  present  in  the  input,  GG  is  not 
read  but  will  be  computed  by  the  solver  itself.  If  the  OPS  card  is  present,  GG  is 
read  and  it  will  be  the  upper  bound  of  the  maximum  eigenvalue  of  the  coefficient 
matrix  using  the  preconditioned  conjugate  gradient  method.  The  default  value  is 
1.0. 


Card  Type 

OPS 

Description 

Maximum  eigenvalue  for  preconditioned  conjugate  gradient  method. 

Repuired 

NO 

Format 

OPS  gg 

Sample 

OPS  1.0 

Field 

Variable 

Value 

Description 

1 

GG 

Upper  bound  for  eigenvalue. 
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Iteration  Parameters 


The  cards  IPl,  JP2,  and  IPS  are  used  to  specify  the  number  of  iterations  for 
the  flow  simulation,  the  transport  simulation,  and  the  coupled  simulation. 


Flow  simulation  (IP1) 

The  IPl  card  is  used  to  specify  the  number  of  iterations  for  solving  the  flow 
equations,  and  the  number  of  cycles  used  to  update  rainfall-seepage  boundary 
conditions,  and  in  determining  convergence  criteria  in  both  steady-state  and 
transient  simulations.  This  card  is  required  for  both  flow-only  and  coupled 
simulations. 


Card  Type 

IPl 

Description 

Iteration  parameters  for  the  flow  simulation. 

Repuired 

YES 

Format 

IP1  niterf  ncylf  npiterf  tolaf  tolbf 

Sample 

IP1  40  10  400  0.0001  0.0001 

Field 

Variable 

Value 

Description 

1 

NITERF 

+ 

Number  of  iterations  allowed  for  solving  the  nonlinear  flow 
equation. 

2 

NCYLF 

+ 

Number  of  cycles  permitted  for  iterating  rainfall  seepage 
boundary  condition  per  time-step. 

3 

NPITERF 

+ 

Number  of  iterations  allowed  for  solving  linearized  flow 
equations  by  pointwise  iterative  solver. 

4 

TOLAF 

+ 

Steady  state  convergence  criterion  for  flow  simulation  FLl. 

5 

TOLBF 

+ 

Transient  convergence  criterion  for  flow  simulation  fLl. 

The  following  iteration  parameters  must  be  designated  for  the  flow 
simulation. 

a.  Number  of  iterations  for  solving  the  nonlinear  flow  equation,  NTTERF. 
Normally  a  value  of  40  is  necessary  for  solving  the  nonlinear  flow 
equations.  However,  if  this  number  is  exceeded,  a  warning  message  is 
issued. 

b.  Number  of  iterations  used  per  time-step  to  check  if  rainfall-seepage 
boundary  conditions  are  properly  used  and  converged,  NCYLF.  If  no 
rainfall-seepage  bound2uy  conditions  are  used,  the  value  should  be  1. 
When  they  are  used,  a  vaJue  of  10  should  be  adequate.  If  10  is  not 
adequate  for  convergence,  a  warning  message  is  issued. 

c.  Number  of  iterations  allowed  for  solving  the  linearized  flow  equation  by 
pointwise  iterative  solver,  NPITERF.  A  value  of  400  is  sufficient  for 
most  problems.  If  this  number  is  exceeded  and  the  solution  does  not 
converge,  the  program  issues  a  warning  message. 
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d.  Steady-state  convergence  criterion  for  the  flow  simulation,  TOLAF. 

This  is  the  absolute  error  allowed  to  determine  if  a  steady-state  solution 
has  converged  for  hydraulic  heads.  A  value  of  0.00001  for  the 
maximum  disturbance  is  usually  sufficient  for  most  problems.  However, 
the  convergent  criterion  should  be  determined  in  the  calibration  of  model 
parameters. 

e.  Transient  convergence  criterion  for  the  flow  simulation,  TOLBF.  This 
is  the  absolute  error  allowed  to  determine  if  hydraulic  heads  have 
converged.  A  value  of  0.0001  for  the  maximum  disturbance  is  usually 
sufficient  for  most  problems.  However,  the  convergent  criterion  should 
be  determined  in  the  calibration  of  model  parameters. 


Transport  simulation  (IP2) 

The  IP2  card  is  used  to  specify  the  number  of  iterations  for  solving  the 
transport  equation,  the  transport  equation  by  pointwise  solver,  and  convergence 
criteria  for  transient  transport  simulations.  This  card  is  required  for  both 
transport  only  and  coupled  simulations. 


Card  Type 

IP2 

Description 

Iteration  parameters  for  the  transport  simulation. 

Required 

YES 

Format 

IP2  nitert  npiterttolbt 

Sampie 

IP2  40  10  400  0.001 

Field 

Variabie 

Vaiue 

Description 

1 

NITERT 

+ 

Number  of  iterations  allowed  for  solving  the  nonlinear  transport 
equation. 

2 

NPITERT 

+ 

Number  of  iterations  allowed  for  solving  linearized  transport 
equation. 

3 

TOLBT 

+ 

Transient  converqence  criterion  for  transport  simulation 

The  following  iteration  parameters  must  be  designated  for  transport 
simulations. 

a.  Number  of  iterations  dlowed  for  solving  the  nonlinear  transport 
equation,  NTTERT.  Normally,  a  value  of  40  is  sufficient.  If  this  number 
is  exceeded  and  the  solution  does  not  converge,  a  warning  message  is 
displayed. 

b.  Number  of  iterations  allowed  for  solving  the  linearized  transport 
equation,  NPITERT.  Normally,  a  value  of  400  is  sufficient  for  solving 
the  linearized  transport  equation  with  the  successive  point  iterative 
solver.  Exceeding  this  number  will  indicate  a  nonconvergent  solution 
and  cause  a  message  to  be  displayed.  When  this  occurs,  a  larger  value 
should  be  used. 
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c.  Transient  convergence  criterion  for  transport  simulation,  TOLBT.  This 
is  the  relative  error  for  determining  if  concentrations  have  converged 
during  transient  simulations.  A  value  of  0.001  is  recommended. 


Coupled  simulation  (IP3) 

The  IP3  card  is  used  to  specify  number  of  iterations  used  in  solving  the 
coupled  flow  and  transport  equations.  This  card  is  used  in  addition  to  the  IP  1 
and  IP2  cards. 


Card  Type 

IPS 

Description 

Iteration  parameters  for  the  coupled  simulation. 

Repaired 

YES 

Format 

IP3  nitfit  omeftt  epss  epst 

Sample 

IPS  10  0.5  0.01  0.05 

Field 

Variable 

Value 

Description 

1 

NITFIT 

+ 

Number  of  iterations  allowed  for  solving  the  coupled  nonlinear 
equations  for  transient  solutions. 

2 

OMEFTT 

+ 

Iteration  parameter  for  solving  the  coupled  nonlinear  equations 
for  transient  solutions. 

3 

EPSS 

+ 

Convergence  criterion  for  head  for  solving  the  coupled 
nonlinear  equations  for  transient  solutions  fLl. 

4 

EPST 

+ 

Convergence  criterion  for  concentration  for  solving  the  coupled 
nonlinear  equations  for  transient  solutions  fM/L^. 

The  following  iteration  parameters  must  be  designated  for  coupled 
simulation. 

a.  Number  of  iterations  allowed  for  solving  the  coupled  nonlinear 
equations  for  transient  solutions,  NTTFTT.  Normally,  a  value  of  10 
should  be  sufficient.  If  this  number  is  exceeded  and  the  solution  does 
not  converge,  a  warning  message  will  be  issued. 

b.  Iteration  parameter  for  solving  the  coupled  nonlinear  equations  for 
transient  solutions,  OMEFTT.  This  parameter  is  the  weighting  factor 
used  with  the  present  and  previous  results  for  solving  the  coupled 
nonlinear  equation  of  the  transient  solution.  A  value  of  0.5  should  be 
used  for  most  problems. 

c.  Convergence  criterion  for  head  for  solving  the  coupled  nonlinear 
equations  for  transient  solutions,  EPSS.  This  is  the  absolute  error 
allowed  for  determining  if  a  coupled  flow  and  transport  solution  for 
hydraulic  head  has  converged.  A  value  of  0.01  for  Ae  maxirrmm 
disturbance  should  be  sufficient. 

d.  Convergence  criterion  for  concentration  for  solving  the  coupled 
nonlinear  equations  for  transient  solutions,  EPST.  This  is  the  relative 
error  allowed  for  determining  if  a  coupled  flow  and  transport  solution 
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for  concentration  has  converged.  A  value  of  0.05  for  the  maximum 
disturbance  should  be  sufficient. 


Particle  Tracking  Parameters 

Particle  tracking,  as  its  name  implies,  is  a  means  of  using  numerical  results 
to  track  fictitious  individual  particles  across  a  numerical  model  mesh.  In  order  to 
accurately  track  particles  over  large  elements  with  large  velocity  gradients,  it 
sometimes  necessary  to  subdivide  the  individual  elements  into  smaller 
subelements. 


Card  Type 

PT1 

Description 

Particle  trackinq  parameters. 

Repuired 

YES 

Format 

PT1  nxw  nvw  nzw  idetq 

Sample 

PT1 1  1  1  1 

Field 

Variable 

Value 

Description 

1 

NXW 

+ 

The  number  of  qrids  for  element  trackinq  in  the  x-direction. 

2 

NYW 

+ 

The  number  of  qrids  for  element  trackinq  in  the  v-direction. 

3 

NZW 

+ 

The  number  of  qrids  for  element  trackinq  in  the  z-direction. 

4 

IDETQ 

Index  of  particle  tracking  pattern. 

1 

Average  velocity  is  used. 

2 

Sinqie  velocity  of  the  startinq  point  is  used. 

The  particle  tracking  parameters  are  entered  on  the  PTl  card.  The  following 
parameters  must  be  specified: 

a.  The  number  of  grids  for  element  tracking  in  the  x-direction,  NXW.  This 
parameter  specifies  how  many  subelements  are  needed  in  the  x-direction. 
The  higher  the  velocity  variation  in  the  x-direction,  the  more 
subelements  are  needed. 

b.  The  number  of  grids  for  element  tracking  in  the  y-direction,  NYW.  This 
parameter  specifies  how  many  subelements  are  needed  in  the  y-direction. 
The  higher  the  velocity  variation  in  the  y-direction,  the  more 
subelements  are  needed. 

c.  The  number  of  grids  for  element  tracking  in  the  z-direction,  NZW.  This 
parameter  specifies  how  many  subelements  are  needed  in  the  z-direction. 
The  higher  the  velocity  variation  in  the  z-direction,  the  more 
subelements  are  needed. 

d.  The  particle  tracking  pattern,  IDETQ.  Two  options  are  available: 

(1)  Average  velocity  is  used  (IDETQ=1).  The  use  of  average  velocity  is 
more  accurate,  and  it  requires  fewer  subelements. 
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(2)  Single  velocity  of  the  starting  point  is  used  (IDETQ=2).  This  option 
should  be  used  when  the  velocity  pattern  is  so  complicated  that  the 
use  of  the  average  velocity  would  fail  to  locate  a  fictitious  particle. 

It  should  be  used  when  a  quick  tracking  is  needed. 


Time  Control  Parameters 


The  TCI  and  TC2  cards  are  used  to  specify  the  total  simulation  time  and  the 
time-step  interval. 


Maximum  simulation  time  (TCI) 

This  is  actual  length  of  time  to  be  simulated.  Once  the  simulation  time 
reaches  the  maximum  simulation  time,  the  simulation  will  be  terminated. 


Card  Type 

TCI 

Description 

Maximum  simuiation  time. 

Repaired 

YES 

Format 

TCI  tmax 

Sample 

TC1  1000 

Field 

Variable 

Value 

Description 

1 

TMAX 

+ 

The  maximum  simulation  time  m. 

Time-step  interval  (TC2) 


The  computational  time-step  can  be  specified  either  as  a  constant  value  or  a 
series  of  time-step  sizes.  The  variable  time-step  should  be  used  when  the 
boundary  conditions  are  changing  rapidly. 


Card  Type 

TC2  ~ 

Description 

Time-step  size. 

Repaired 

YES 

Format 

TC2  idt  delMdtxv 

Sample 

TC2  0  1.0 

Field 

Variable 

Value 

Description 

1 

IDT 

0 

Time-step  type. 

Constant. 

1 

Variable. 

2 

DELT 

+ 

If  IDT=0.  constant  time-step. 

IDTXY 

+ 

If  IDT=1 ,  index  of  xy  card  for  variable  time-step  values.  The 
x  values  in  the  series  represent  times  at  which  the  time-step 
size  will  change.  The  y  values  represent  the  time-step 
sizes. 
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The  XY  Series  format  (XY1) 


The  TC2  card  in  the  previous  section  uses  the  XY  Series  format  (XYl  card) 
to  specify  how  the  time-step  size  changes  with  time.  The  XYl  card  is  used  by 
many  of  the  FEMWATER  input  cards  as  a  consistent  means  of  designating  a  list 
of  pairs  of  numbers.  This  can  be  thought  of  as  the  x-  and  y-coordinates  of  a 
curve,  although  the  values  do  not  necessarily  have  to  correspond  to  a  curve.  In 
many  cases,  the  x- value  represents  a  time  and  the  y-value  represents  some 
parameter  that  is  changing  with  time.  However,  the  series  can  represent  any 
relationship.  Each  XYl  card  has  an  index  that  is  referenced  by  other  cards. 

XYl  cards  can  appear  anywhere  in  the  input  file  and  are  often  listed  at  the  end  of 
the  file. 


Card  Type 

XY1 

Description 

XY  Series.  Used  to  define  a  sequence  of  pairs  of  numbers. 

Repaired 

YES 

Format 

XY1  i  n  dx  dy  rep  begc  tname 

XI  yi 

X2y2 

XnVn 

Sample 

XYl  1  5  0  0  0  0  head 

4.0  0.0 

4.1  2.0 

4.2  7.0 

4.3  8.0 

4.4  9.5 

Fieid 

Variabie 

Vaiue 

Description 

1 

i 

+ 

index  of  xv  series. 

2 

N 

+ 

The  number  of  xy  pairs  in  the  series. 

3 

DX 

Vaiue  is  ignored  by  FEMWATER. 

4 

DY 

Vaiue  is  ignored  by  FEMWATER. 

5 

REP 

Vaiue  is  ignored  bv  FEMWATER. 

6 

BEGC 

Vaiue  is  ignored  by  FEMWATER. 

7 

TNAME 

text 

A  character  string  representing  the  name  of  the  XY  series. 

8 

X,Y 

±,  ± 

After  the  XYl  card,  the  xy  values  of  the  series  are  listed, 
one  pair  per  line,  up  to  N  pairs. 

Output  Control  Parameters 

There  are  two  basic  types  of  output  generated  by  FEMWATER.  The  first 
type  consists  of  printed  text  output  summarizing  the  input  data,  the  progress  of 
the  simulation  (convergence  criteria,  etc.),  and  a  summary  of  the  results.  This 
type  of  output  is  printed  to  the  screen  and  to  the  printed  output  file  described  in 
Chapter  2.  The  other  type  of  output  consists  of  a  series  of  binary  or  text  solution 
files  that  can  be  used  for  graphical  post-processing.  The  format  of  the  solution 
files  is  described  in  Appendix  C.  The  output  control  card  group  (OCl,  OC2, 
OC3,  and  OC4)  is  used  to  specify  what  information  is  to  be  written  to  the  printed 
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output  file  and  saved  to  the  solution  files.  The  cards  also  control  the  interval  at 
which  the  output  is  printed  or  saved. 


Print  interval  (OC1) 

The  OCl  card  controls  the  frequency  at  which  data  are  written  to  the  printed 
output  file  and  whether  or  not  diagnostic  output  and  the  cyclic  change  of  rainfall 
seepage  nodes  are  printed. 


Card  Type 

OC1 

Description 

interval  for  writing  to  printed  output  file. 

Required 

YES 

Format 

OCl  ibug  ichng  iopt  kprt/nprint 

Sample 

OC1  0  0  01  /*  Print  at  every  time-step  */ 

OCl  0  015  /*  Print  at  specified  times  ’/ 

1.0 

2.0 

4.0 

8.0 

16.0 

Field 

Variable 

Value 

Description 

1 

IBUG 

0 

1 

Do  not  print  diagnostic  output. 

Print  diagnostic  output. 

2 

ICHNG 

0 

1 

Do  not  print  cyclic  change  of  rainfall/seep,  nodes. 

Print  cyclic  change  of  rainfall/seep,  nodes. 

3 

JOPT 

0 

1 

Print  at  specified  interval. 

Print  at  specified  set  of  time  values. 

4 

KPRT 

KPRINT 

+ 

+ 

Print  interval  if  JOPT=0,  or 

Total  number  of  specified  times  if  JOPT=1 . 

5 

+ 

If  JOPT=1 ,  after  the  OCl  card,  list  the  specified  print 
times,  one  per  line,  up  to  KPRINT  times. 

Dia^ostic  output  (BBUG).  The  integer,  IBUG,  is  an  indicator  for 
diagnostic  output  control.  The  diagnostic  output  will  help  users  locate  errors  in 
the  input  data. 

Cyclic  change  of  rainfall-seepage  nodes  (ICHNG).  The  integer,  ICHNG, 
is  an  indicator  for  controlling  the  printout  of  the  convergence  behavior  at 
rainfall-seepage  nodes. 

Print  frequency  (JOPT).  The  information  that  is  written  to  the  printed 
output  file  can  either  be  written  out  at  a  regular  interval  or  at  a  specified  set  of 
time  values.  The  time  values  in  the  specified  set  of  time  values  should 
correspond  to  computational  time-steps. 
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Print  options  (OC2) 

The  OC2  card  is  used  to  specify  what  information  is  written  to  the  printed 
output  file.  The  information  is  written  at  the  interval  specified  in  the  OCl  card. 


Card  Type 

OC2 

Description 

Options  for  writing  to  printed  output  file. 

Repuired 

YES 

Format 

002  nselt  koro(i) 

Sample 

0C2  314  7 

Field 

Variable 

Value 

Description 

1 

NSELT 

+ 

Total  number  of  print  options  selected.  This  field  should 

by  followed  by  NSELT  values,  each  of  which  ranges  from 

0-7  as  explained  below. 

2+ 

KPRO(I) 

0 

Print  nothing. 

1 

Print  flow  information  at  the  boundary. 

2 

Print  total  head. 

3 

Print  pressure  head. 

4 

Print  concentration. 

5 

Print  flux. 

6 

Print  nodal  moisture  content. 

7 

Print  Darcy  velocity. 

The  following  items  can  be  printed: 

a.  Flow  information  at  the  boundary.  The  rate  of  change,  incremental,  and 
total  flow  through  the  boundary  are  printed. 

b.  Total  head.  The  total  head  at  each  node  is  printed. 

c.  Pressure  head.  The  pressure  head  at  each  node  is  printed. 

d.  Concentration.  The  concentration  at  each  node  is  printed. 

e.  Flux.  The  flux  at  each  node  is  printed. 

f.  Moisture  content.  The  moisture  content  at  each  node  is  printed. 

g.  Darcy  velocity.  The  Darcy  velocity  is  printed  at  each  node. 

Save  interval  (OC3) 

The  OC3  card  controls  the  frequency  at  which  the  computed  solution  is 
saved  to  output  files  for  post-processing.  The  solution  can  be  saved  at  a  regular 
interval  or  at  a  specified  set  of  time  values.  The  time  values  in  the  specified  set 
of  time  values  should  correspond  to  computational  time-steps. 

The  OC3  is  also  used  to  specify  whether  the  files  should  be  saved  in  text  or 
binary  format.  The  binary  format  results  in  much  smaller  file  sizes,  an  issue 
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which  may  be  very  important  when  ranning  transient  analyses  on  large  meshes. 
The  text  format  requires  more  disk  space,  but  it  can  be  viewed  with  a  text  editor 
and  transferred  between  different  computer  platforms.  Binary  files,  however, 
cannot  be  transferred  between  different  computer  platforms. 


Card  Type 

"oci 

Description 

Interval  for  savinq  solution  files  for  post-processinq. 

Required 

YES 

Format 

OC3  ifile  kopt  kdsk/npost 

Sampie 

OC3  0  0  1  /’  Save  at  every  time-step  */ 

OC3  0  15  r  Save  at  specified  times  */ 

1.0 

2.0 

4.0 

8.0 

16.0 

Fieid 

Variable 

Value 

Description 

1 

IFILE 

0 

1 

Save  in  text  format. 

Save  in  binary  format. 

2 

KOPT 

0 

1 

Save  at  specified  interval. 

Save  at  specified  set  of  time  values. 

3 

KDSK 

KPOST 

+ 

+ 

Save  interval  if  KOPT=0,  or 

Total  number  of  specified  times  if  KOPT=1. 

4+ 

+ 

If  KOPT=1 ,  after  the  OC3  card,  list  the  specified  save 
times,  one  per  line,  up  to  KPOST  times. 

Save  options  (OC4) 

The  OC4  card  is  used  to  specify  what  information  is  saved  to  the  solution 
files  for  post-processing  or  for  use  as  initial  conditions  for  subsequent 
simulations.  Each  set  of  information  is  written  to  a  separate  file,  which  is  either 
binary  or  text  depending  on  the  option  selected  in  the  OC3  card.  The  formats  of 
the  solution  files  are  described  in  Appendix  C.  The  solutions  are  saved  one 
time-step  at  a  time  at  the  interval  specified  in  the  OC3  card.  At  each  time-step, 
one  value  is  written  for  each  node. 
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Card  Type 

OC4 

Description 

Ootions  for  savinq  solution  files  for  oost-processing. 

Reauired 

YES 

Format 

OC4  kselt  ksave(i) 

Sample 

004  3  13  4 

Field 

Variable 

Value 

Description 

1 

KSELT 

+ 

Total  number  of  solution  files  selected.  This  field  should 

be  followed  by  KSELT  values,  each  of  which  ranges  from 

0-5  as  explained  below. 

2+ 

KSAVE(I) 

0 

Save  nothing. 

1 

Save  pressure  head. 

3 

Save  nodal  moisture  content. 

4 

Save  velocity. 

5 

Save  concentration. 

The  following  types  of  solution  files  can  be  saved; 

a.  Pressure  head. 

b.  Nodal  moisture  content  computed  at  nodes  as  an  average  of  surrounding 
elements. 

c.  Darcy  velocity. 

d.  Concentration. 
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Introduction 

The  second  of  the  four  primary  groups  of  information  in  the  model  file  is  the 
material  properties.  The  material  properties  include  fluid  properties  and  soil 
properties.  One  set  of  fluid  properties  is  defined  for  the  entire  mesh,  and  one  set 
of  soil  properties  is  defined  for  each  of  the  soil  or  aquifer  types  referenced  by  the 
element  material  ID’s. 


File  Format 

The  set  of  cards  in  the  model  file  corresponding  to  material  properties  is 
shown  in  Figure  5. 


MPl  kcp 

/♦  Cond.  Or  perm.  Flag*/ 

MP2  i  kxx  kxy  kz  kxy  kxz  kyz  alpha  por 

/*  Cond.  Or  perm.  Values  •/ 

MP3  iho  vise  grav  betap 

/*  Dens.  And  vise.  Of  water  */ 

MP4  al  a2  a3  a4  a5  a6  a?  a8 

/*  Dens.  And  vise.  Coeff.s  */ 

MP5  I  k  gamma  al  at  am  t  decay  n  deckw  decks 

/*  Disp./diffusion  coeff.s  */ 

SPl  I  ihm  ihe  ihw 

/*  Soil  prop.For  unsat.  zone  */ 

Figure  5.  The  material  properties  cards  in  the  model  file 


Fluid  Properties 

The  fluid  property  cards  input  to  the  model  file  are  used  to  specify  the 
density,  the  viscosity  and  the  compressibility  of  fluid,  and  the  acceleration  of 
gravity.  The  acceleration  of  gravity  is  not  a  fluid  property,  but  it  is  input  on  the 
same  line  as  the  density  and  viscosity  for  convenience.  Specifying  the 
acceleration  of  gravity  allows  the  user  to  use  any  desired  units  for  the  fluid 
properties,  soil  properties,  and  all  other  parameters  input  to  FEMWATER.  All 
parameters  should  be  consistent  with  the  units  of  the  specified  acceleration  of 
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gravity.  For  example,  if  the  gravity  constant  is  specified  in  units  of  m/day  ,  all 
length  units,  including  mesh  coordinates,  must  be  in  meters  and  all  time  units, 
including  time-step  data,  must  be  specified  in  units  of  days.  However,  since  the 
most  common  mass  concentration  unit  used  in  groundwater  is  milligrams  per 
liter  (mg/^  )  in  SI  units,  the  concentrations  in  FEMWATER  should  be  specified 
in  units  of  mg/^  (ppm). 

FEMWATER  can  be  used  to  model  density-driven  flow  and  transport.  Thus, 
relationships  must  be  defined  between  concentration,  density,  and  viscosity.  The 
relationships  used  by  FEMWATER  are: 

=  +  (31) 

Po 

and 

i=a,+a,C+a,C=+a.C’  (32) 

where 

Po,  Ho  =  density  and  viscosity  of  fresh  water 

a,...ag  =  parameters  used  to  define  concentration  dependence  of  water 
density  and  viscosity 

C  =  chemical  concentration 

Thus,  values  of  Po,  Ho,  and  Uj . .  .ag  must  be  specified  by  the  user. 


Density  and  viscosity  of  fresh  water  (MP3) 

The  density,  viscosity,  and  compressibility  of  fresh  water ,  po.  Ho,  P’>  atid  the 
acceleration  of  gravity  are  specified  with  the  MP3  card. 


Card  Type 

MP3 

Description 

Density  and  viscosity  of  fresh  water. 

Repaired 

YES 

Format 

MP3  rho  vise  qrav  betap 

Sample 

MP3  1000.0  4.68  1.3e+8  0.0 

Field 

Variable 

Value 

Description 

1 

RHO 

+ 

Density  of  fresh  water  fM/L’l. 

2 

Vise 

+ 

Dynamic  viscosity  of  water  fM/L/n. 

3 

GRAY 

+ 

Acceleration  of  gravity  fL/n. 

4 

BETAP 

+ 

Comoressibility  of  water  fl/Ll. 
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Concentration  dependence  coefficients  (MP4) 

The  coefficients,  ai.-.ag,  which  are  used  to  define  the  concentration 
dependence  of  water  density  and  viscosity,  are  specified  with  the  MP4  card.  For 
density-independent  flow,  the  parameters  ai  and  as  are  set  to  1.0,  and  the 
parameters  a2,  as,  a4,  a^,  arj,  and  ag  are  set  to  zero. 

For  density-dependent  flow,  some  or  all  of  the  parameters  a2-a4  and  ag-ag 
should  be  nonzero.  For  an  example,  if  the  density  depends  on  concentration 
only,  the  parameter  ai  is  set  to  1.0,  a2is  set  to  0.001,  ag,  84  are  set  to  zero,  and  the 
unit  of  concentration  is  kg/Iiter  (ppt). 


Card  Type 

MP4 

Description 

Concentration  dependence  coefficients. 

Repaired 

YES 

Format 

MP4  a1  a2  a3  a4  aS  a6  a7  a8 

Sampie 

MP4  1 .0  0.0  0.0  0.0  1 .0  0.0  0.0  0.0 

Fieid 

Variable 

Value 

Descriptior) 

1 

RHOMUd) 

1 

Coefficient  at . 

2 

RHOMU(2) 

± 

Coefficient  a2.  The  units  are  the  inverse  of  concentration 
units. 

3 

RHOMU(3) 

± 

Coefficient  a3.  The  units  are  the  inverse  of  concentration 
units  squared. 

4 

RHOMU{4) 

± 

Coefficient  a4.  The  units  are  the  inverse  of  concentration 
units  cubed. 

5 

RHOMU(5) 

1 

Coefficient  a5. 

6 

RHOMU(6) 

± 

Coefficient  a6.  The  units  are  the  inverse  of  concentration 
units  squared. 

7 

RHOMU(7) 

± 

Coefficient  a7.  The  units  are  the  inverse  of  concentration 
units  cubed. 

8 

RHOMU(8) 

± 

Coefficient  a8. 

Soil  Properties 

Each  element  in  the  mesh  is  assigned  a  material  corresponding  to  the  zone  or 
aquifer  in  which  the  element  is  located.  The  material  ID  is  an  index  to  a  list  of 
soil  properties.  The  list  of  soil  properties  is  input  to  the  model  file.  The  soil 
properties  that  must  be  specified  are  hydraulic  conductivity,  compressibility, 
porosity,  dispersion,  diffosion,  decay  coefficients,  the  first  order  biodegradation 
rate  of  dissolved  and  absorbed  phase,  and  three  water  retention  curves  describing 
how  moisture  content,  relative  conductivity,  and  water  capacity  vary  with 
pressure  head  in  the  unsaturated  zone. 
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Hydraulic  conductivity  (MP1,  MP2) 

The  hydraulic  conductivity  for  each  soil  type  is  specified  using  the  MPl  and 
MP2  cards.  The  MPl  card  is  used  to  specify  whether  the  conductivity  values 
should  be  interpreted  as  hydraulic  conductivity  or  as  permeability. 


Card  Type 

MPl 

Description 

Hydraulic  conductivity  versus  oermeability. 

Repuired 

YES 

Format 

MP1  kcp 

Sampie 

MPl  0 

Field 

Variable 

Value 

Description 

1 

KCP 

0 

Hydraulic  conductivity. 

1 

Permeability. 

The  hydraulic  conductivity  K  is  a  symmetric  tensor  of  secondorder.  It 
relates  the  flux  F  to  the  gradient  of  total  head,  H,  in  a  rectangular  system  (x,y,z) 
as 


dH 

dH 

dH 

1 

II 

K 

.K.-) 

dH 

dH 

dH 

dy 

+  Ky.—  ) 

dH 

dH 

dH 

Fz  =  -(K„^ 

dy 

+  Kzz-^  ) 

Since  the  tensor  is  symmetric,  a  coordinate  system  (^,V,0  can  be  found  such  that 
(Long  1961): 

dH 


dH 


F?  =  -K5t 


dH 

dC 


(34) 


Therefore,  the  hydraulic  conductivity  tensor  K  can  be  written  as 


K., 

K„' 

1 _ 

0 

0 

K  = 

= 

0 

0 

K., 

K. 

0 

0 

(35) 


The  transformation  of  K  between  the  coordinate  system  (^,V,Q  and  (x,y,z) 
is  given  as  follows. 
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Let  6,  <p,  and  (p  be  the  Euler  angles  between  the  and  (x,y,z) 

coordinates:  6  is  the  angle  between  the  and  z-axis  on  the  plane  perpendicular 
to  the  y-axis,  (j)  is  the  angle  between  ^  and  x-axis  on  the  plane  perpendicular  to  z- 
axis,  and  9  is  the  third  angle  (Morse  and  Feshbach  1978).  The  relationship 
between  (x,y,z)  and  (4,V,C)  is  given  by 


X 

«1  «2  «3 

y 

>  — 

A  A  A3 

•« 

¥ 

z 

V.  J 

ri  72  73. 

in  which 

tti  =  sin(p  sitKj)  +  cos(p  cos^  cosd 
0.2  =  cos^  sin^  -  sincp  cos<j)  cos6 
0.2  =  cos  <j>  cos  6 

Pi  =  sin  <p  coscj)  -  sin(()  sinip  cosO 
p2  =  cos(p  cos<p  +  sin(p  sin^  cosd 
Ps  =  -sin<p  sin  6 
Yi  =  -cos(psin6 
Y2  =  sin<p  sinO 
Ys  =  cosd 


(36) 


where  a’s,  p’s,  and  y’s  are  the  directional  cosines  between  two  coordinate 
systems.  These  directional  cosines  are  symmetrical  with  respect  to  the  two 
coordinate  systems.  The  transformation  of  the  flux  vector  between  two 
coordinate  systems  is  similar  to  that  between  the  coordinates,  i.e.. 


A 

' 

'Py 

>  — 

AAA 

* 

F,- 

A. 

Jt  77  7t. 

F,. 

Substituting  Equations  33  and  34  into  Equation  37  gives 


(37) 
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K., 

yx 

yy 

K. 

K., 

K 

K 

K, 


yz 


dx 

/y  fy 

dH 

tXi  1X2  tXj 

dy 

*  “ 

A  A  A 

dH 

.7i  72  73. 

.  dz 

0 

0 


0 

K 

w 

0 


0 

0 


K, 


dH 


dH 


dxjf 

dH 


di; 


(38) 


By  the  chain  rule. 


dH 

dx  dy  dz 

dH 

^  ^  ^ 

dx 

dH 

dx  dy  dz 

dH 

dy/ 

dy/  dy/  dy/ 

dH 

dx  dy  dz 

dH 

.  dz . 

The  Jacobian  can  be  obtained  from  Equation  37: 


dx  dy  dz 

^  d^  d^ 

dx  dy 
dyr  dy/  dy/ 
dx  dy  dz 
^  ^  ^ 


a, 

A  A  A 

71  72  73. 


(39) 


(40) 


Substituting  Equation  40  into  Equation  39,  then  the  resulting  equation  into 
Equation  38  gives 


1 _ 

«! 

1 

is 

Kyy 

Ky. 

= 

A 

A 

A 

.7i 

72 

73. 

L 

0 

r 

0 


0 

0 


^2 


A 

A 

A 


71 

72 

73 


(41) 


Relabeling  the  directional  cosines  in  Equation  41  gives 
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721 

0 

0 

7i2 

Kyy 

= 

712 

722 

732 

0 

0 

721 

722 

723 

7i3 

723 

733. 

0 

0 

K.. 

.731 

732 

733. 

in  which 


(42) 


Yii  =  sin(psin(l)  +  cos(pcos<pcosd 

Yi2  =  sin(j>cos^  -  cos(p  sin(j)  cos9 

Yi3  =  -cos(psin6 

Y21  =  cos(psin<j>  -  sin(pcos<j>cosd 

Y22  =  cos^cos^  +  sin(psin<pcosd 

Y23  =  sincpsin  6 

Y31  =  sinOcos^ 

Y32  =  -sin6sin(j) 

Y33  =  cos6 


Finally  the  working  equations  to  compute  the  hydraulic  conductivity  components 
in  the  (x,y,z)  coordinate  are  given  by  the  following  equation,  when  the  principal 
components  in  the  coordinate  and  the  three  Euler  angles  between  the 

(x,y,z)  and  (^,V,Q  coordinates  are  known; 


K„  = 

7i^i 

+ 

721 

731 

Kxy  = 

7ii 

7i2 

+ 

721  722 

K,,ny  4-  731 

732 

Kxz  = 

7ii 

7i3 

+ 

721  723 

Kw  731 

733 

% 

Kyx  = 

7i2 

7ii 

+ 

722  721 

Kyv  4-  732 

731 

Kyy  = 

7i2 

^  + 

722 

K\|njf  4" 

73^2 

Ky.= 

7i2 

7i3 

+ 

722  7a 

Kyi);  4"  732 

733 

K„  = 

7i3 

7ii 

+ 

723  721 

733 

731 

% 

K.y  = 

7i3 

7i2 

% 

+ 

723  722 

733 

732 

K?? 

7i3 

+  ' 

7a 

Kyy  4" 

733 

For  anisotropic  cases,  the  three  principal  components  should  be  all  distinct, 
i.e.,  ^  ^  The  nine  components  of  hydraulic  conductivity  tensor  in 

the  (x,y,z)  coordinate  system  must  be  computed  given  the  three  principal 
components  and  three  Euler  angles. 

a.  For  isotropic  cases,  the  three  principal  components  should  be  all  equal, 
i.e.,  Kxx  =  Kyy  =  Kzz  =  K.  All  off-diagonal  terms  of  the  hydraulic 
conductivity  tensor  must  be  zero  no  matter  what  coordinate  system  is 
used. 
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b.  For  horizontal  isotropic  cases,  all  off-diagonal  terms  must  be  zero,  Kxx 
and  Kyy  must  be  equal  to  K,  and  K^z  must  be  different  from  K. 


c.  For  the  vertical  isotropic  cases,  all  off-diagonal  terms  must  be  zero,  Kyy 
and  Kzz  must  be  equal  to  K,  and  Kxx  must  be  different  from  K. 

The  hydraulic  conductivity  values,  compressibility,  and  porosity  of  the 
medium  are  entered  on  the  MP2  card.  The  values  should  correspond  to  either 
hydrauHc  conductivity  or  permeability  depending  on  the  status  of  the  KCP 
variable  on  the  MPl  card. 


Card  Type 

MP2 

Description 

Hydraulic  conductivity  tensor. 

Required 

YES 

Format 

MP2 1  kxx  kyy  kzz  kxy  kxz  kyz  alphap  por 

Sample 

MPl  10  0.003  0.003  0.003  0.0  0.0  0.0  0.0  0.0 

Field 

Variable 

Value 

Description 

1 

I 

+ 

Material  ID. 

Kxx 

PROPF(1,l) 

+ 

Saturated  xx  hydraulic  conductivity  [L/T],  or  saturated 
peimeabilitv  fL*l  of  the  medium  i. 

Kyy 

PROPF(2.l) 

+ 

Saturated  yy  hydraulic  conductivity  [L/T],  or  saturated 
permeability  [L^  of  the  medium  i. 

Kzz 

PROPF(3,l) 

+ 

Saturated  zz  hydraulic  conductivity  [L/T],  or  saturated 
permeability  [L^  of  the  medium  i. 

Kxy 

PROPF(4,l) 

+ 

Saturated  xy  hydraulic  conductivity  [L/T],  or  saturated 
permeability  fL*1  of  the  medium  i. 

Kxz 

PROPF(5,l) 

+ 

Saturated  xz  hydraulic  conductivity  [L/T],  or  saturated 
permeability  [L*]  of  the  medium  I. 

Kyz 

PROPF(6.l) 

+ 

Saturated  yz  hydraulic  conductivity  [L/T],  or  saturated 
permeability  [L'l  of  the  medium  I. 

ALPHAP 

PROPF(7,l) 

+ 

The  modified  compressibility  of  the  medium  i  [l/L]. 

POR 

PROPF(8,l) 

+ 

The  effective  porosity  of  the  medium  i  [unitlessl. 

Dispersion/diffusion  coefficients  (MP5) 

The  MP5  card  is  used  to  specify  the  bulk  density,  the  parameters  required  to 
define  the  isotherm  relationship,  the  dispersion  coefficients,  the  tortuosity,  the 
decay  constant,  and  the  first-order  biodegradation  rate  constant  through 
dissolved  and  adsorbed  phases  for  each  soil  type. 
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Card  Type 

MP5 

Description 

Dispersion/diffusion  coefficients. 

Repuired 

YES 

Format 

MP5  i  k  aamma  al  at  am  t  decay  n  declw  decks 

Sample 

MP5  1 0  0.0  1 200.0  5.0  5.0  0.0  1 .0  0.0  0.0  0.0  0.0 

Field 

Variable 

Value 

Description 

i 

I 

+ 

Material  type  ID. 

k 

PROPT(1,l) 

+ 

Distribution  coefficient,  K„,  Freundlich  K,  or  Langmuir  K  for 
medium  i  [L’/Ml. 

qamma 

PROPT(2,l) 

+ 

Bulk  density  for  medium  i  fM/L^l. 

al 

PROPT(3,l) 

+ 

Lonqitudinal  dispersion  for  medium  i  fLI. 

at 

PROPT(4,l) 

+ 

Transverse  dispersion  for  medium  i  fLl. 

am 

PROPT{5,l) 

+ 

Molecular  diffusion  coefficient  for  medium  i  FLYTI. 

t 

PROPT(6,l) 

+ 

Tortuosity  for  medium  i  funitlessl. 

decay 

PROPT(7.l) 

+ 

Decay  constant  for  medium  i  ri/r|. 

n 

PROPT(8,l) 

+ 

Freundlich  n  or  Lanqmuir  S_„  for  medium  i  funitlessl. 

deckw 

PROPT(11,l) 

+ 

The  first-order  biodegradation  rate  constant  through 
dissolved  phase  for  medium  i  fl/TI. 

decks 

PROPT(12,l) 

+ 

The  first-order  biodegradation  rate  constant  through 
adsorbed  phase  for  medium  i  fl/Tl. 

The  density  is  specified  as  a  bulk  density,  in  units  of  [M/L^].  The  bulk 
density  of  a  soil  corresponds  to  its  oven-dried  mass  divided  by  its  in  situ  volume. 

The  isotherm  is  used  to  define  a  relationship  between  the  amount  of 
contaminant  adsorbed  by  the  soil,  S  (mass  of  contaminant  /  mass  of  soil)  and  the 
concentration  of  the  contaminant,  C.  Three  types  of  isotherms  can  be  used  in 
FEMWATER: 

a.  Linear  isotherm: 


S  =  K,C 


(44) 


b.  Langmuir  isotherm: 

.  S^KC 
1-hKC 


(45) 


c.  Freundlich  isotherm: 


S  =  Kr 


(46) 


where 

Kd  =  distribution  coefficient  [L^/M] 
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Smax  =  maximum  concentration  of  medium  in  the  Langmuir  nonlinear 
isotherm 

K  =  coefficient  in  the  Langmuir  or  Freundlich  isotherm 

n  =  power  index  in  the  Freundlich  nonlinear  isotherm. 

Note  that  a  linear  isotherm  can  be  obtained  using  Equation  44  by  setting  n  equal 
to  1.  The  type  of  isotherm  to  be  used  is  specified  with  the  KSORP  variable  on 
the  OPS  card.  However,  one  set  of  coefficients,  Ka,  Smax.  K,  and  n,  is  specified 
for  each  material  type  on  the  MP5  card. 

Diffusion  and  dispersion  are  modeled  in  FEMWATER  using  the  following 
relationship: 


W 

0D  =  a, I Vl6  +  (a^  -  a, )  -j^  +  a, 0x6 


(47) 


where 

6  =  moisture  content 

D  =  dispersion  coefficient  tensor 

ot  =  transverse  dispersion  [L] 

d  =  Kronecker  delta  tensor 

IVI  =  the  magnitude  of  the  Darcy  velocity  V  [L/T] 

Cl  =  longitudinal  dispersion  [L] 

Om  =  molecular  diffusion  coefficient  [L^/T] 

T=  tortuosity 

Thus,  the  parameters,  ot,  ai,  a„,  and  t  must  be  specified  by  the  user.  They  can  be 
obtained  by  experiment  or  literature  review. 

A  decay  constant,  X  [1/T],  must  also  be  specified  for  each  medium.  The 
decay  constant  is  used  to  simulate  biodegradation  or  radioactive  decay  of  the 
contaminant. 


Soil  properties  for  unsaturated  zone  (SPI) 

The  final  set  of  parameters  that  must  be  specified  for  each  medium  is  a 
sequence  of  three  curves  defining  how  the  moisture  content,  relative 
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conductivity,  and  water  capacity  vary  as  a  function  of  pressure  head  in  the 
unsaturated  zone. 


Card  Type 

SP1 

Description 

Soil  properties  for  unsaturated  zone. 

Repaired 

YES 

Format 

SP1  i  ihm,  ihc,  ihw 

Sample 

SP1  10  345 

Field 

Variable 

Value 

Description 

1 

1 

+ 

Material  type  index. 

2 

IHM(I) 

+ 

Index  of  moisture  content  versus  pressure  head  XY  series. 

3 

IHC(I) 

+ 

Index  of  relative  conductivity  versus  pressure  head  XY  series. 

4 

IHW(I) 

+ 

Index  of  water  capacity  versus  pressure  head  XY  series. 

Relative  conductivity.  The  governing  equation  for  flow  of  water  through  a 
variably  saturated  porous  medium  has  effective  hydraulic  conductivity  and 
storage  terms.  The  effective  hydraulic  conductivity  can  be  rewritten  as  the 
product  of  nonlinear  and  constant  terms  in  the  form: 


K(h)  =  K,Ks  (48) 

where  Kr  is  the  relative  conductivity,  ranging  in  value  from  0.0  to  1.0,  and  Ks  is 
the  saturated  hydraulic  conductivity.  The  change  of  relative  hydraulic 
conductivity  is  caused  by  changes  in  moisture  content,  resulting  in  the 
preferential  movement  of  water  through  certain  pathways,  due  to  the  influence  of 
capillary  forces.  As  the  soil  becomes  less  saturated,  the  flow  of  water  becomes 
restricted  to  the  pore  sequences  of  smaller  radii.  This  results  in  a  reduction  in 
the  spatially  averaged  effective  hydraulic  conductivity. 

Moisture  content.  Moisture  content  in  the  unsaturated  zone  is  a  function  of 
the  pressure  head.  The  more  negative  the  pressure  head,  the  lower  the  moisture 
content.  The  curve  defining  the  relationship  between  moisture  content  and 
pressure  head  should  vary  between  the  saturated  moisture  content,  6s,  and  the 
residual  moisture  content,  6r.  The  saturated  moisture  content  is  equal  to  the 
porosity  of  the  medium  since  all  of  the  void  space  is  filled  with  fluid.  Under 
unsaturated  conditions,  however,  some  of  the  void  space  is  filled  with  air;  thus, 
the  moisture  content  is  less  than  the  medium’s  porosity.  The  residual  moisture 
content  represents  the  amount  of  water  that  cannot  be  removed  from  a  soil  by 
gravity  drainage  (even  under  large  suction  pressure)  because  it  adheres  to  the 
grains  of  the  soil. 

Water  capacity.  The  water  capacity  curve  is  equal  to  the  slope  of  the 
moisture  content  versus  pressure  head  curve.  Although  this  curve  could  be 
determined  automatically  by  FEMWATER  from  the  moisture  content  curve,  it  is 
input  by  the  user  to  avoid  errors  resulting  from  the  approximate  piecewise  linear 
nature  of  the  moisture  content  curve. 
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Generating  the  curves.  Ideally,  the  relative  conductivity,  moisture  content, 
and  water  capacity  curves  are  determined  directly  by  performing  a  series  of  tests 
on  the  soils  involved  in  the  study.  However,  in  many  cases  they  can  be 
approximated  using  a  set  of  measured  or  approximated  constants  and  a  set  of 
empirical  relationships.  For  example,  one  option  for  generating  the  curves  is  to 
use  the  van  Genuchten  functions  (van  Genuchten  1980): 

K,  =ef[l-(l-0f)"]'  (49) 


and 


ft  = 


i+W)' 


-y 


for  h  <  0 


(50) 


0,  =1  forh>0 


where 


ew=e.+0e(e,-e.) 


(51) 


(52) 


Y  = 


(53) 


and 


0^  =  moisture  content  (dimensionless) 

0g  =  effective  moisture  content  (dimensionless) 

05=  saturation  moisture  content  (dimensionless) 

05  =  residual  moisture  content  (dimensionless) 
j3,Y  =  soil-specific  exponents  (dimensionless) 
a=  soil-specific  coefficient  (1/L) 

Table  3  lists  a  set  of  saturated  and  residual  moisture  contents  and  the  van 
Genuchten  a  and  j3  terms  for  a  variety  of  common  soil  types.  When  applying  the 
a  term,  care  should  be  taken  to  convert  it  to  the  proper  units. 
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Table  3 

Representative  Soil  Parameters 

Soil  Type 

Saturated 
Moisture 
Content  0. 

Residual 
Moisture 
Content  0, 

a[cm'] 

Clav' 

0.38 

0.068 

0.008 

1.09 

Clav  loam 

0.41 

0.095 

0.019 

1.31 

Loam 

0.43 

0.078 

0.036 

1.56 

Loam  sand 

0.41 

0.057 

0.124 

2.28 

Silt 

0.46 

0.034 

0.106 

1.37 

Silt  loam 

0.45 

0.067 

0.020 

1.41 

Siltv  clav 

0.36 

0.070 

0.005 

1.09 

Siltv  clav  loam 

0.43 

0.089 

0.010 

1.23 

Sand 

0.43 

0.045 

0.145 

2.68 

Sandv  clav 

0.38 

0.100 

0.027 

1.23 

Sandv  clav  Loam 

0.39 

0.100 

0.059 

1.48 

Sandv  loam 

0.41 

0.065 

0.075 

1.89 

Note:  *  Agricultural  soil,  less  than  60%  clay 

Source:  Carsel  and  Parrish  (1988) 

Model  convergence.  The  unsaturated  zone  soil  property  curves  can  have  a 
significant  effect  on  model  convergence.  If  any  portion  of  a  curve  has  a  sharp 
change  in  gradient  (extreme  curvature),  FEMWATER  may  have  a  difficult  time 
converging  during  the  iterative  solution  process.  If  convergence  is  a  problem, 
careful  editing  of  the  curves  can  often  resolve  the  problem. 

Entering  the  curves.  The  user  must  supply  the  moisture  content,  relative 
conductivity,  and  water  capacity  versus,  pressure  head  curves  in  tabular  form. 
The  curves  are  specified  at  the  end  of  the  model  file  using  a  set  of  xy  series 
cards.  The  curves  are  referred  to  on  the  SPl  card  by  the  ID’s  of  the  xy  series. 
The  format  of  the  xy  series  card  is  described  in  Chapter  4,  “The  XY  Series 
format.” 
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6  Boundary  Conditions 


introduction 

The  boundary  conditions  supported  by  FEMWATER  include  point  sources 
and  sinks^ ,  Dirichlet  (nodal),  Cauchy  (flux),  Neumann  (flux  gradient),  and 
variable  boundaries  (rainfall/evaporation  and  seepage).  In  each  case,  the 
prescribed  values  can  be  either  constant  or  vary  with  time.  Not  assigning 
boundary  conditions  to  a  mesh  boundary  is  equivalent  to  a  no-flux  boundary  in 
flow  simulation  and  to  a  no-diffusion  flux  boundary  in  transport  simulation. 


Choosing  Appropriate  Boundary  Conditions 

FEMWATER  offers  the  user  several  types  of  boundary  conditions  in  an 
attempt  to  cover  the  range  of  valid  boundary  conditions  that  can  be  encountered 
in  finite  element  groundwater  modeling.  Some  confusion  over  which  boundary 
conditions  are  appropriate  for  various  situations  will  imdoubtedly  occur.  The 
following  discussion  is  intended  to  help  in  choosing  appropriate  boundary 
conditions  in  FEMWATER  simulations. 

The  most  common  boundary  condition  used  in  groundwater  modeling  is  the 
specified  head  (Dirichlet)  boundary  condition.  If  heads  are  known  at  a  given 
location  over  the  length  of  the  simulation,  this  is  the  most  appropriate  boundary 
condition  to  use. 

As  a  general  rale,  the  variable  boundary  condition  should  be  used  when 
modeling  rainfall  and  evaporation.  This  boundary  condition  is  the  most  robust 
because  it  allows  FEMWATER  the  freedom  to  change  the  boundary  condition  to 
a  specified  head  boundary  condition  for  both  over-  and  undersaturated 
conditions  (see  the  discussion  of  this  boundary  condition  later  in  this  chapter). 
This  feature  is  especially  applicable  in  transient  simulations  where  the  water 
table  may  fluctuate  over  the  simulation  time  span. 


'  Point  sources/sinks  are  technically  not  boundary  conditions.  However,  they  are  included  in  this 
chapter  for  convenience. 
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Some  modelers  may  be  more  accustomed  to  using  the  specified  flux 
(Cauchy)  boundary  condition  for  infiltration/evaporation  and  inflow/outflow 
flux  conditions.  Indeed,  in  some  instances  it  is  possible  to  know  the  flux  on 
element  faces  with  a  reasonable  degree  of  assurance.  However,  the  user  is 
cautioned  that  the  specified  flux  is  always  taken  in  a  direction  parallel  to  the 
normal  of  the  element  face  to  which  this  boundary  condition  is  applied  (see  the 
discussion  of  this  boundary  condition  later  in  this  chapter).  Therefore,  in  a 
simulation  where  a  known  flux  condition  exists,  the  specified  flux  value  must  be 
adjusted  for  all  element  faces  whose  normals  deviate  from  the  known  flux 
direction.  In  most  groundwater  modeling  scenarios,  and  particularly  in  transient 
simulations,  the  model  is  best  suited  either  by  using  specified  head  boundary 
conditions  (for  inflow/outflow  conditions  on  model  boundaries)  or  by  allowing 
the  numerical  model  to  change  a  flux  boundary  condition  to  a  constant  head 
condition  if  conditions  warrant  (for  a  infiltration/evaporation  flux  condition). 

The  variable  boundary  condition  in  FEMWATER  provides  the  user  with  the 
latter  capability. 

Gradient  flux  (Neumann)  boundary  conditions,  as  indicated  in  the  discussion 
of  this  boundary  condition  later  in  this  chapter,  are  rarely  encountered  in  natural 
conditions  but  are  included  in  FEMWATER  for  completeness  and  for  research 
applications. 


File  Format 

Boundary  conditions  are  stored  in  the  model  file.  The  set  of  model  file  cards 
corresponding  to  boundary  conditions  is  shown  in  Figure  6. 


•PSl  nodeid  flowseries 

/*  Point  source,  flow  rate  */ 

PS2  nodeid  concseries 

/*  Point  source,  concentration  ♦/ 

DBl  nodeid  headseries 

/*  Dirichlet,  head  */ 

DB2  nodeid  concseries 

/*  Dirichlet,  concentration  */ 

CBl  elemid  faceid  fluxseries 

/*  Cauchy,  flux  rate  */ 

CB2  elemid  faceid  concseries 

/*  Cauchy,  chemical  flux  rate  */ 

RSI  elemid  faceid  fluxseries 

/*  Variable,  flux  */ 

RS2  elemid  faceid  concseries 

/*  Variable,  concentration  */ 

RS3  hcon  hmin 

/♦  Variable,  pond  depth,  min  hd.  */ 

Figure  6.  The  boundary  condition  cards  in  the  model  file 


Element  Faces 


The  Cauchy,  Neumann,  and  variable  boundary  conditions  described  in  later 
sections  are  all  “flux”  type  boundary  conditions  and  are  assigned  to  element 
faces.  An  element  face  is  referenced  by  two  indices,  the  element  number  and  a 
face  index.  The  element  number  is  simply  the  element  ID  The  numbering 
scheme  used  to  determine  face  indices  is  shown  in  Figure  7. 
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Figure  7.  Numbering  scheme  for  element  faces 


Point  Sources/Sinks  (PS) 

Point  source/sink  type  boundary  conditions  consist  of  a  flow  rate  assigned  to 
nodes  in  the  mesh.  Sources  and  siidts  typically  correspond  to  injection  or 
extraction  wells.  If  the  well  is  an  injection  well,  the  concentration  of  the  fluid 
being  injected  can  also  be  specified.  The  flow  rate  is  a  positive  value  for  an 
injection  well  and  is  a  negative  value  for  an  extraction  well. 

Point  sources/sinks  are  specified  with  the  PSl  and  PS2  cards.  The  PSl  card 
is  used  to  specify  the  flow  rate  (L^/T),  and  the  PS2  card  is  used  to  specify  the 
concentration  (W  L^).  Both  cards  contain  a  field  which  is  an  index  to  an  XY 
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Series  card.  The  XY  Series  card  (specified  elsewhere  in  the  file)  describes  how 
the  values  (flow  rate  and  concentration)  vary  with  time.  If  the  value  is  constant 
(does  not  vary  with  time),  only  one  value  should  be  given  in  the  XY  Series.  The 
XY  Series  card  is  described  in  Chapter  4. 


Card  Type 

PS1 

Description 

Point  source/sink  flow  rate  for  both  flow  and  transport  simulation. 

Repaired 

NO 

Format 

PS1  nodeid  fiowseries 

Sample 

PS1  324  1 

Field 

Variable 

Value 

Description 

nodeid 

NPWF(I) 

+ 

Node  number. 

fiowseries 

IWTYPF(I) 

+ 

Index  of  the  flow  rate  XY  Series. 

Card  Type 

PS2 

Description 

Point  source/sink  concentration  for  transport  simulation. 

Repaired 

NO 

Format 

PS2  nodeid  concseries 

Sample 

PS2  324  2 

Field 

Variable 

Value 

Description 

nodeid 

NPWT(I) 

+ 

Node  number. 

concseries 

IvmPT(l) 

+ 

Index  of  the  concentration  XY  Series. 

Dirichlet  Boundary  Conditions  (DB) 

A  Dirichlet  boundary  condition  consists  of  a  prescribed  total  head  at  a  node. 
A  concentration  may  also  be  specified.  The  Dirichlet  boundary  condition  is 
usually  applied  to  soil-water  interface  such  as  stream,  artificial  impoundments, 
lakes,  and  coastal  lines. 

The  governing  flow  equations  for  FEMWATER  are  written  in  terms  of 
pressure  head.  Thus,  the  formal  definition  of  a  Dirichlet  boundary  condition  for 
head  is: 


h  =  h/Xb,yb,Zb,t)  onBj, 


(54) 


where  (X(,,yj,,z^)  is  the  spatial  coordinate  on  the  boundary,  is  the  Dirichlet 
boundary,  and  h^  is  the  prescribed  pressure  head.  Since  it  is  more  convenient  to 
specify  total  head  rather  than  pressure  head,  Dirichlet  boundary  conditions  are 
input  as  total  head  values,  and  they  are  converted  internally  by  FEMWATER  to 
pressure  heads  by  subtracting  the  elevation  of  the  nodes  where  the  boundary 
conditions  are  assigned. 

Thus,  the  formal  definition  of  a  Dirichlet  concentration  boundary  condition 
for  head  is: 
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C  =  Cd(Xb,yb,Zb)  onBd, 


(55) 


where  Cd  is  the  prescribed  concentration  at  the  Dirichlet  boundary. 

Dirichlet  boundary  conditions  are  specified  with  the  DB 1  and  DB2  cards. 
The  total  head  (L)  is  specified  with  the  DB  1  card  for  flow  simulation,  and  the 
concentration  (M/L^)  is  specified  with  the  DB2  card  for  transport  simulation.  In 
both  cases,  the  values  are  specified  with  XY  Series  cards.  If  the  head  or 
concentration  does  not  change  with  time,  only  one  value  should  be  listed  in  the 
XY  Series. 


Card  Type 

DB1 

Description 

Dirichlet  boundary  conditions,  prescribed  total  head  for  flow  simulation. 

Reauired 

NO 

Format 

DB1  nodeid  headseries 

Sample 

DB1  324  2 

Field 

Variable 

Value 

Description 

nodeid 

NPDBF(I) 

+ 

Node  number. 

headseries 

IDTYPF(I) 

+ 

Index  of  the  head  XY  Series. 

Card  Type 

DB2 

Description 

Dirichlet  boundary  conditions,  prescribed  concentration  for  transport  simulation. 

Reouired 

NO 

Format 

DB2  nodeid  concseries 

Sample 

DB2  324  3 

Field 

Variable 

Value 

Description 

nodeid 

NPDBTd) 

+ 

Node  number. 

concseries 

JDTYPTd) 

+ 

Index  of  the  concentration  XY  Series. 

Cauchy  Boundary  Conditions  (CB) 


A  Cauchy  boundary  condition  consists  of  a  fluid  flux  prescribed  at  a 
boundary  element  face.  The  concentration  of  the  fluid  may  also  be  specified. 
Formally,  the  Cauchy  boundary  condition  for  flow  can  be  stated  as: 


n  K 


^Vh  +  Vz 


V 


=  qc(Xb.yb’Zb’t)  onB,, 


J 


(56) 


where 

n  =  outward  unit  vector  normal  to  the  boundary 
(Xb,  yb,  Zb)  =  spatial  coordinate  on  the  boundary 
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K  =  hydraulic  conductivity  tensor 

po  =  reference  (clean)  fluid  density 

p  =  solution  density 

h  =  pressure  head 

z  =  elevation 

qc  =  flux  rate 

Be  =  Cauchy  boundary 

The  Cauchy  boundary  for  transport  can  be  stated  as; 

n-(vC-0D-Vc)  =  q,(xj,,yb,Zb,t)  onB^,  (57) 


where 

V  =  Darcy  velocity 
C  =  concentration 
B  =  moisture  content 
D  =  dispersion  coefficient  tensor 
qc  =  prescribed  flux  rate 
Be  =  Cauchy  boundary 

The  Cauchy  boundary  condition  is  typically  applied  to  surface  water  bodies 
with  known  infiltration  rates  through  the  layers  of  the  bottom  sediments  or  liners 
into  the  subsurface  media.  The  flux  rate  for  flow  and  transport  is  a  negative 
value  when  it  leaves  the  system  and  is  a  positive  value  when  it  enters  the  system. 

Cauchy  boundary  conditions  are  specified  with  the  CBl  and  CB2  cards.  The 
flux  rate  (L/T)  is  specified  with  the  CBl  card  for  flow  simulation,  and  the 
material  flux  rate  (mass  per  unit  area  per  time,  M/L^/T)  is  specified  with  the  CB2 
card  for  transport  simulation.  In  both  cases,  the  values  are  specified  with  XY 
Series  cards.  If  the  flux  rate  does  not  change  with  time,  only  one  value  should  be 
listed  in  the  XY  Series. 
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Card  Type 

CB1 

Description 

Cauchv  boundary  conditions,  prescribed  flux  rate  for  flow  simulation. 

Repaired 

NO 

Format 

CB1  elemid  faceid  fluxseries 

Sampie 

CB1  353  3  8 

Fieid 

Variabie 

Vaiue 

Description 

elemid 

NPCBF(I) 

+ 

Element  number. 

faceid 

IDCF(I) 

+ 

Element  face  ID. 

fluxseries 

ICTYPF(I) 

+ 

Index  of  the  flux  XY  Series. 

Card  Type 

CB2 

Description 

Cauchv  boundary  conditions,  prescribed  flux  rate  for  transport  simulation. 

Repuired 

NO 

Format 

CB2  elemid  faceid  concseries 

Sampie 

CB2  353  3  9 

Fieid 

Variable 

Vaiue 

Description 

elemid 

NPCBT(I) 

+ 

Element  number. 

faceid 

IDCT(I) 

+ 

Element  face  ID. 

concseries 

ICTYPT(I) 

+ 

Index  of  the  concentration  XY  Series. 

Neumann  Boundary  Conditions  (NB) 

The  Neumann  boundary  condition  rarely  occurs  in  nature,  but  is  included  for 
completeness.  It  consists  of  a  fluid  flux  gradient  prescribed  at  a  boundary 
element  face.  The  concentration  of  the  fluid  may  also  be  specified.  Formally, 
the  Neumann  boundary  condition  for  flow  can  be  stated  as: 

-n-K-y  Vh  =  q„(Xb,yb,Zb,t)  onB„,  (58) 

where 

n  =  outward  imit  vector  normal  to  the  boundary 

(Xb,  yb,  Zb)  =  spatial  coordinate  on  the  boundary 

K  =  hydraulic  conductivity  tensor 

Po  =  reference  (clean)  fluid  density 

p  =  solution  density 

h  =  pressure  head 

Qn  =  prescribed  flux  rate 
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Bn  =  Neumann  boundary 


The  Neumann  boundary  for  transport  can  be  stated  as: 

n-(-eDVC)  =  q„(xb,yb,Zb,t)  onB„,  (59) 


where 

C  =  concentration 
9  =  moisture  content 
D  =  dispersion  coefficient  tensor 
qn  =  prescribed  flux  rate 
Bn  =  Neumaim  boundary 

The  flux  gradient  for  flow  and  transport  is  a  negative  value  when  it  leaves 
the  system  and  is  a  positive  value  when  it  enters  the  system. 

Neumann  boundary  conditions  are  specified  with  the  NBl  and  NB2  cards. 
The  flux  gradient  (L/T)  is  specified  with  the  NB 1  card  for  flow  simulation,  and 
the  flux  rate  (mass  per  unit  area  per  time,  M/L^/T)  is  specified  with  the  NB2  card 
for  transport  simulation.  In  both  cases,  the  values  are  specified  with  XY  Series 
cards.  If  the  flux  gradient  does  not  change  with  time,  only  one  value  should  be 
listed  in  the  XY  Series. 


Card  Type 

NB1 

Description 

Neumann  boundary  conditions,  prescribed  flux  gradient  for  flow  simulation. 

Repuired 

NO 

Format 

NB1  elemid  faceid  fluxseries 

Sample 

NBl  353  3  8 

Field 

Variable 

Value 

Description 

elemid 

NPNBF(I) 

+ 

Element  number. 

faceid 

iDNF(l) 

+ 

Element  face  ID. 

fluxseries 

INTYPF(I) 

Index  of  the  flux  gradient  XY  Series. 

Card  Type 

NB2 

Description 

Neumann  boundary  conditions,  prescribed  flux  rate  for  transport  simulation. 

Required 

NO 

Format 

NB2  elemid  faceid  fluxseries 

Sample 

NB2  353  3  9 

Field 

Variable 

Value 

Description 

elemid 

NPNBT(I) 

+ 

Element  number. 

faceid 

IDNT(I) 

+ 

Element  face  ID. 

fluxseries 

INTYPT(I) 

+ 

Index  of  the  concentration  XY  Series. 

56 


Chapter  6  Boundary  Conditions 


Variable  Boundary  Conditions  (RS) 


Variable  boundary  conditions  are  typically  applied  to  the  top  faces  of  the 
uppermost  layer  of  elements  in  the  mesh  (i.e.,  at  the  soil-air  interface),  and  are 
used  to  simulate  evaporation  and  seepage  due  to  precipitation.  Variable 
boundary  conditions  are  called  “variable”  not  because  they  can  vary  with  time, 
but  because  they  correspond  to  either  a  Dirichlet  or  a  Cauchy  boundary 
condition  depending  on  the  potential  evaporation,  the  conductivity  of  the  media, 
and  the  availability  of  water  such  as  rainfall. 

During  the  precipitation  period,  variable  boundary  conditions  are  defined  as 
h  =  hp(Xb,yb,Zb,t)  onB^,  (60) 


or 


-n  K 


^Vh  +  Vz 
VP  J 


=  ^p(Xi,,yi,,Zi„t)  onB^, 


(61) 


During  the  nonprecipitation  period,  variable  boundary  conditions  are  defined 


as 


h  =  h  (Xb,yb,Zb,t)  onB,, 


(62) 


or 


h  =  h„(Xb,yb.Zb’t)  onB,, 


(63) 


or 


-n  K 


^Vh-hVz 

Ip  ) 


=  qe(Xb’yb’Zb’t) 


onB,, 


(64) 


where 

h  =  pressure  head 
hp  =  ponding  depth 

(Xb,yb,Zb)  =  spatial  coordinate  on  the  boundary 
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Bv  =  variable  boundary 

n  =  outward  unit  vector  normal  to  the  boundary 

K  =  hydraulic  conductivity  tensor 

Po  =  reference  (clean)  fluid  density 

p  =  solution  density 

z  =  elevation 

qp  =  precipitation 

hm  =  minimum  pressure  head 

qe  =  maximum  evaporation  rate 

Only  one  of  Equations  60-64  is  used  at  any  given  time  on  the  variable 
boundary.  During  the  precipitation  period,  rainfall  infiltrates  to  the  subsurface  at 
a  rate  equal  to  the  prescribed  precipitation  flux,  qp,  according  to  Equation  61.  If 
the  precipitation  flux  exceeds  the  infiltration  capacity  of  the  soil,  the  pressure 
head  is  not  allowed  to  rise  above  the  ponding  depth,  hp.  If  it  does,  the  boundary 
condition  switches  to  the  Dirichlet  condition  of  Equation  60,  i.e.,  the  pressure 
head  is  set  equal  to  the  ponding  depth. 

During  the  nonprecipitation  period,  the  boundary  condition  corresponds  to 
one  of  Equations  62-64.  Once  again,  the  pressure  head  is  not  allowed  to  rise 
above  the  ponding  depth  (Equation  62).  If  the  pressure  head  is  below  the 
ponding  depth,  the  boundary  condition  is  switched  to  the  Cauchy  condition  of 
Equation  64  and  the  flux  is  set  equal  to  the  specified  evaporation  flux,  qe.  If  the 
pressure  head  drops  below  the  minimum  pressure  head,  hn„  the  boundary 
condition  is  switched  to  the  Dirichlet  condition  of  Equation  63,  and  the  pressure 
head  is  set  equal  to  the  minimum  pressure. 

During  periods  when  the  variable  boundary  condition  is  set  to  a  Cauchy 
boimdary  condition,  the  specified  flux  value  is  taken  to  be  directed  vertically  (i.e. 
in  either  the  positive  or  negative  z-direction)  in  order  to  simulate  the  primary 
direction  of  rainfall  and  evaporation.  To  calculate  the  volumetric  flux  entering 
the  system,  FEMWATER  takes  the  specified  flux  rate  and  multiplies  this  value 
by  the  area  of  the  element  face  projected  onto  a  horizontal  xy-plane 
(perpendicular  to  the  flux  direction  ).  For  element  faces  that  deviate  from  the 
horizontal  by  some  degree,  this  calculation  will  ensure  that  the  flux  is  not 
artificially  exaggerated.  The  variable  boundary  condition  differs  from  the 
specified  flux  boundary  condition  in  this  matter.  The  specified  flux  boundary 
condition  always  assumes  that  the  specified  flux  is  pardlel  to  the  element  face 
normal  and  thus  the  normal  element  face  area  is  used  for  the  volumetric  flux 
calculation.  This  feature  of  the  variable  boundary  condition  allows 
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FEMWATER  to  more  accurately  simulate  real  world  conditions  of  rainfall 
infiltration  and  evaporation  on  a  topographically  varying  surface. 

For  transport  simulations,  during  both  the  precipitation  and  nonprecipitation 
periods,  the  variable  boundary  condition  is  either  a  Cauchy  boundary  condition 
with  a  given  total  flux,  or  a  Neumann  with  zero  gradient  flux.  Mathematically, 
this  can  be  expressed  as: 

ii-(vC-eD-VC)  =  n-VC,(xb,yb,Zb,t)  ifn-V<0  (65) 

n-(-0D  VC)  =  O  ifii  V>0  (66) 


where 

C  =  concentration 

V  =  Darcy  velocity  tensor 

6  =  moisture  content 

D  =  dispersion  coefficient  tensor 

Cv  =  concentration  at  the  variable  boundary 

If  the  flow  direction  through  the  boundary  is  into  the  region,  the  Cauchy 
condition  of  Equation  65  governs  and  the  concentration  in  the  incoming  fluid  is 
equal  to  the  prescribed  concentration  Cy.  If  the  flow  direction  through  the 
boundary  is  out  of  the  region,  the  Neumann  condition  of  Equation  66  governs, 
meaning  that  the  concentration  of  the  outgoing  fluid  is  whatever  it  was  before  it 
left  the  region. 

Variable  boundary  conditions  are  defined  by  specifying  the  minimum 
pressure  head,  hm,  the  ponding  depth,  hp,  the  flux  rates  qp  and  qe,  and  the 
concentration,  Cy.  The  minimum  pressure  head  and  the  ponding  depth  are 
specified  on  the  RS3  card.  Only  one  RS3  card  is  entered  in  the  model  file.  The 
specified  minimum  pressure  head  and  ponding  depth  are  used  for  all  element 
faces  where  a  variable  boundary  condition  is  specified. 

The  flux  rates  (L/T)  for  precipitation  and  evaporation,  qp  and  q^  are  both 
specified  on  the  RS 1  card  as  an  index  to  a  single  XY  Series.  The  values  in  the 
series  that  are  positive  are  treated  as  precipitation  flux,  and  the  values  that  are 
negative  are  interpreted  as  evaporation  flux.  Thus,  cycles  of  rainfall/evaporation 
can  be  designated  on  one  curve  by  varying  the  curve  in  the  XY  Series  from 
positive  to  negative  and  vice  versa. 
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The  concentration  (M/L^)  is  entered  on  the  RS2  card  as  an  index  to  an  XY 
series.  The  ponding  depth  (L)  and  the  minimum  pressure  head  (L)  are  both 
specified  on  the  RS3  card. 


Card  Type 

RS1 

Description 

Variable  boundary  condition,  rainfall/evaporation/seepage  flux  rate  for  flow 
simulation. 

Reauired 

NO 

Format 

RS1  elemid  faceid  fluxseries 

Sample 

RS1  8749  4  12 

Field 

Variable 

Value 

Description 

elemid 

NPVBF(I) 

+ 

Element  number. 

faceid 

IDRF(I) 

+ 

Element  face  index. 

fluxseries 

IVTYPF(I) 

Index  number  of  rainfall/evaporation  time  series  for  the 
element  face. 

Card  Type 

RS2 

Description 

Variable  boundary  condition,  rainfall/seepaqe  concentration  for  transport  simulation. 

Repuired 

NO 

Format 

RS2  elemid  faceid  concseries 

Sample 

RS2  8749  413 

Field 

Variable 

Value 

Description 

elemid 

NPVBT(I) 

+ 

Element  number. 

faceid 

IDRT(I) 

Element  face  index. 

concseries 

IVTYPT(I) 

+ 

Index  number  of  concentration  time  series  for  the 
element  face. 

Card  Type 

RS3 

Description 

Variable  boundary  condition,  ponding  depth  and  minimum  pressure  head. 

Repuired 

NO 

Format 

RS3  hcon  hmin 

Sample 

RS3  0.35  -8.0 

Field 

Variable 

Value 

Description 

1 

HCON 

+ 

The  ponding  depth. 

2 

HMIN 

± 

The  minimum  pressure  head. 

60 


Chapter  6  Boundary  Conditions 


7  Initial  Conditions 


Introduction 

Whenever  a  FEMWATER  analysis  is  performed,  a  set  of  initial  conditions 
must  be  defined.  Initial  conditions  define  the  initial  status  of  the  pressure  head 
and  concentration.  Which  initial  conditions  are  required  for  a  particular  problem 
depends  on  the  type  of  simulation  that  is  being  performed.  The  rules,  options, 
and  formats  for  defining  initial  conditions  for  FEMWATER  are  described  in  this 
chapter. 


Types  of  Initial  Conditions 

Three  types  of  initial  conditions  are  possible  for  a  FEMWATER  simulation: 
cold  starts,  hot  starts,  and  flow  solutions.  Cold  starts  are  used  to  establish  a  set 
of  initial  values  at  the  begiiming  of  a  steady-state  or  transient  simulation.  Hot 
starts  are  used  to  continue  a  previous  run  of  FEMWATER  without  having  to 
start  over  from  the  beginning.  Flow  solutions  are  used  to  define  the  flow  field 
that  is  necessary  when  performing  a  transport-only  simulation  (as  opposed  to 
coupled  flow  and  transport). 


Cold  Starts 

FEMWATER  uses  an  iterative  method  to  find  a  solution  to  flow  and 
transport  problems.  An  iterative  approach  involves  starting  with  an  initial  set  of 
heads  and  concentrations  and  utilizing  an  iterative  solver  to  gradually  modify  the 
initial  values  until  they  converge  to  a  set  of  values  that  satisfy  the  underlying 
governing  equations.  This  initial  set  of  values  is  called  a  cold  start. 


Steady  state  versus  transient 
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Cold  starts  must  be  specified  regardless  of  whether  the  simulation  is  steady- 
state  or  transient.  However,  the  approach  taken  to  define  the  cold  start  may  vary 
depending  on  whether  the  simulation  is  steady  state  or  transient.  With  a  steady- 
state  simulation,  the  set  of  initial  values  chosen  is  often  not  critical.  If  the  values 
are  closer  to  the  final  solution,  the  simulation  will  converge  with  fewer  iterations 
(resulting  in  a  much  shorter  execution  time).  However,  most  sets  of  initial 
conditions  will  converge  to  the  same  final  solution. 

With  a  transient  simulation,  the  specification  of  the  initial  condition  becomes 
more  critical.  In  this  case,  the  initial  condition  should  represent  the  actual  state 
of  the  aquifer  at  the  begiiming  of  the  simulation  period.  In  order  for  the 
simulation  to  be  accurate,  the  initial  condition  must  be  compatible  with  the 
specified  boundary  conditions  and  stresses.  For  example,  suppose  that  the 
transient  effect  of  a  new  injection/extraction  weU  system  on  an  aquifer  is  to  be 
simulated.  If  an  incompatible  initial  condition  is  used,  the  response  of  the 
aquifer  in  the  early  part  of  the  simulation  will  be  due  in  part  to  the  wells  and  in 
part  to  an  adjustment  of  the  aquifer  that  is  required  to  make  the  heads  and 
concentrations  compatible  with  the  other  stresses  and  boundary  conditions. 
Therefore,  it  is  often  best  to  first  perform  a  steady  state  simulation  with  the 
boundary  conditions  that  will  not  change  during  the  transient  simulation,  and 
then  use  the  solution  to  this  simulation  as  the  initial  condition  for  the  transient 
simulation.  For  this  example,  this  would  involve  first  running  a  steady-state 
simulation  with  all  of  the  boundary  conditions  except  for  the  injection  and 
extraction  wells.  The  solution  to  Ais  simulation  would  then  be  used  as  the  initial 
condition  for  the  transient  simulation,  and  the  response  of  the  aquifer  would  be 
due  solely  to  the  added  stresses  from  the  injection  and  extraction  wells. 

In  the  case  of  concentration,  a  steady-state  simulation  is  often  not  applicable 
since  transport  seldom  achieves  steady-state  conditions.  In  many  cases,  the 
initial  condition  wiU  correspond  to  an  existing  plume  that  has  been  characterized 
via  three-dimensional  interpolation  from  measurements  at  scattered  locations. 


Required  values 

If  a  flow-only  simulation  is  performed,  a  set  of  pressure  heads  is  required  for 
the  cold  start  initial  condition.  If  a  transport-only  simulation  is  performed,  a  set 
of  concentrations  is  required  (in  addition  to  the  flow  solution  as  explained 
below).  If  a  coupled  flow  and  transport  simulation  is  being  performed,  a  set  of 
heads  and  a  set  of  concentrations  are  both  required. 


Convergence 

The  values  chosen  for  a  cold  start  can  have  a  significant  impact  on 
convergence,  ff  the  selected  cold  start  differs  greatly  from  the  final  solution,  i.e., 
is  highly  incompatible  with  the  model  stresses,  FEMWATER  may  have 
difficulty  converging.  The  quickest  and  most  robust  solution  is  achieved  when 
cold  start  values  are  as  close  to  the  final  solution  as  possible. 
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Hot  Starts 


Since  FEMWATER  is  a  nonlinear,  finite  element  model,  solutions  can 
require  extensive  computation  time  even  for  moderately  sized  meshes, 
particularly  when  a  transient  simulation  involving  numerous  time-steps  is  being 
performed.  During  a  modeling  study,  it  is  not  uncommon  to  encounter  situations 
where  the  solution  progresses  well  to  a  certain  point  in  time  and  then  diverges  or 
“blows  up.”  This  simation  can  often  be  solved  by  decreasing  the  time-step  size 
near  the  point  in  time  where  the  solution  degraded.  In  such  cases,  it  is  often 
useful  to  be  able  to  restart  the  simulation  at  some  point  just  before  the  problem 
occurred  rather  than  starting  over  from  the  begmning.  In  order  to  do  this,  the 
“hot  start”  type  of  initial  condition  is  required. 

Another  application  of  hot  start  fiOies  is  “spinning  up”  a  model.  As 
mentioned  in  the  previous  section,  if  the  cold  start  initial  condition  differs 
greatly  from  the  final  solution,  the  stresses  due  to  the  boundary  conditions  may 
cause  FEMWATER  to  oscillate  and  have  difficulty  converging.  One  solution  to 
this  problem  is  to  select  a  better  cold  start.  Another  solution  is  to  start  with  a  set 
of  boundary  conditions  that  are  more  compatible  with  the  specified  cold  start  and 
compute  a  solution.  The  boundary  conditions  are  then  modified  slightly  and  the 
solution  from  the  first  simulation  is  used  as  a  hot  start  for  the  second  simulation. 
This  process  is  repeated  and  the  boundary  conditions  are  gradually  changed  until 
they  correspond  to  the  desired  values. 


Hot  start  time 

A  hot  start  consists  of  a  set  of  files  representing  a  transient  solution  of 
pressure  head,  moisture  content  (optional),  velocity  (optional),  and  concentration 
and  a  hot  start  time.  FEMWATER  reads  the  transient  solutions  and  locates  the 
set  of  values  at  the  time-step  matching  the  hot  start  time.  This  set  of  values  is 
then  loaded  into  FEMWATER  as  initial  conditions  and  the  simulation  is 
restarted. 

When  the  initial  or  original  simulation  is  performed,  the  time-variant 
boundary  conditions  are  typically  defined  on  a  time  scale  whose  beginning 
corresponds  to  the  start  time  of  the  simulation.  When  a  hot  start  is  performed,  it 
is  not  necessary  to  redefine  time- variant  boundary  conditions  so  that  they  are 
defined  from  a  point  corresponding  to  the  beginning  of  the  new  start  time. 
FEMWATER  uses  the  specified  hot  start  time  and  automatically  interpolates  the 
time-variant  boundary  conditions  to  the  correct  beginning  time. 


Required  values 

The  solution  files  necessary  for  a  hot  start  depend  on  the  type  of  simulation. 
If  a  flow-only  simulation  is  being  performed,  pressure  head  is  required.  If  a 
transport-only  simulation  is  being  performed,  concentration  and  pressure  head 
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are  required.  If  a  coupled  flow  and  transport  simulation  is  being  performed, 
pressure  head  and  concentration  are  required.  The  velocity  and  nodal  moisture 
content  are  used  for  post-processing  only. 


File  format 

The  files  used  for  hot  starts  correspond  to  the  data  set  solution  files  output 
by  FEMWATER  for  post-processing.  These  files  consist  of  solutions  organized 
by  time-step  and  are  sequential  lists  of  scalar  or  vector  values,  with  one  value 
listed  per  node.  The  nodal  values  are  used  for  post-processing.  The  output 
control  options  which  control  data  set  output  are  described  in  the  section  “Save 
interval  (OC3).”  The  formats  of  the  data  set  files  are  described  in  Appendix  C. 


Flow  Solutions 

A  third  type  of  initial  condition  is  required  when  a  transport-only  simulation 
is  being  performed.  A  transport-only  simulation  uses  a  previously  computed 
flow  solution  (steady-state  or  transient)  to  define  the  three-dimensional  flow 
field  required  to  properly  model  the  contaminant  migration.  The  flow  solution 
consists  of  two  options;  pressure  head  or  velocity  and  pressure  head  written  to 
data  set  files.  When  the  first  option  is  used,  the  needed  velocity  and  moisture 
content  for  transport  simulation  are  computed  from  the  pressure  head.  When  the 
second  option  is  used,  the  needed  moisture  content  is  computed  from  the 
pressure  head.  Thus,  when  miming  the  flow  simulation,  care  should  be  taken  to 
ensure  that  the  proper  options  are  selected  for  solution  output  (see  “Save  options 
(OC4)”). 

The  flow  solution  for  a  transport-only  simulation  is  used  in  combination  with 
either  a  cold  start  or  a  hot  start.  With  a  cold  start,  a  set  of  initial  concentration 
values  is  provided  for  concentration  in  addition  to  the  steady-state  or  transient 
flow  solution.  With  a  hot  start,  a  transient  concentration  solution  and  a  hot  start 
time  is  provided  in  addition  to  the  flow  solution. 


Summary  of  Initial  Condition  Requirements 

It  can  sometimes  be  confusing  to  determine  which  combination  of  initial 
conditions  is  required  for  a  particular  situation.  The  sets  of  initial  conditions 
required  for  all  possible  scenarios  are  shown  in  Table  4.  The  X’s  in  the  lower 
right  comer  of  the  table  indicate  which  sets  of  values  are  required.  For  example, 
for  a  cold  start  with  a  coupled  flow  and  transport  simulation,  pressure  head  and 
concentration  must  be  specified. 
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Table  4 

1  Initial  Condition  Requirements 

Initial  Condition  I 

Simulation  Type 

Press.  Head 

Moisture  Cent. 

Velocity 

Concentration  | 

1  Cold  Start  I 

Flow 

X 

Transport 

X 

X 

Coupled 

X 

X 

I  Hot  Start  [ 

Flow 

X 

X' 

X’ 

X 

Transport 

X 

X’ 

X 

Coupled 

X 

X’ 

X’ 

X 

1  '  Solution  files  are  optional  1 

Initial  conditions  are  specified  in  a  two-step  process.  First  of  all,  a  set  of 
cards  must  be  entered  in  the  model  file  specifying  the  basic  control  parameters 
for  initial  conditions  such  as  the  hot  start  time  and  the  format  of  the  initial 
condition  files  (text  or  binary).  Next,  the  initial  conditions  are  assembled  in  a  set 
of  data  set  files.  The  data  set  files  should  be  located  in  the  same  directory  as  the 
other  input  files.  A  set  of  cards  are  entered  in  the  super  file  that  specify  the 
names  of  the  data  set  files  containing  the  initial  conditions. 


Model  File  input 

A  set  of  input  cards  is  required  in  the  model  file  to  designate  the  control 
parameters  for  the  initial  conditions.  The  cards  are  summarized  in  Figure  8. 
Each  of  the  cards  is  explained  in  more  detail  in  the  following  sections. 


ICS  istart 

/*  Cold  start  vs.  hot  start  */ 

ICT  hstime 

/*  Hot  start  time  */ 

ICC  icon  conval 

r  Use  constant  concentration  */ 

ICH  ihead  hconst 

r  Compute  press  from  const  head  */ 

ICF  icfile  ivfile  iffu 

/*  Initial  cond.  file  format  */. 

Figure  8.  The  initial  condition  cards  in  the  model  file 


Start  type  (ICS) 


The  ICS  card  is  used  to  specify  whether  the  simulation  is  to  be  started  with  a 
cold  start  or  a  hot  start. 
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Card  Type 

ICS 

Description 

Initial  condition  start  type. 

Required 

YES 

Format 

ICS  istart 

Sample 

ICS  1 

Field 

Variable 

Value 

Description 

1 

ISTART 

0 

Cold  start. 

1 

Hot  start. 

Hot  start  time  (ICT) 

If  a  hot  start  is  specified  on  the  ICS  card,  the  hot  start  time  must  be  specified  on 
the  ICT  card.  This  time  is  used  to  interpolate  the  time  variant  boundary 
conditions  to  a  proper  starting  point  and  to  locate  the  proper  set  of  initial 
conditions  from  the  multiple  time-steps  defined  in  the  data  set  files  of  pressure 
head,  moisture  content,  etc. 


Card  Type 

ICT 

Description 

Initial  condition  hot  start  time. 

Required 

YES 

Format 

ICT  hstime 

Sample 

ICT  234.0 

Field 

Variable 

Value 

Description 

1 

HSTIME 

± 

Hot  start  time. 

Constant  or  variable  concentration  (ICC) 

When  defining  a  single  set  of  concentration  values  for  a  cold  start  initial 
condition,  it  is  often  useful  to  use  a  constant  value  of  concentration  everywhere 
in  the  problem  domain.  For  example,  in  many  cases,  an  initial  condition  of  zero 
concentration  everywhere  in  the  problem  domain  is  appropriate.  The  ICC  card 
can  be  used  to  easily  define  a  constant  concentration  for  the  entire  mesh.  This 
alleviates  the  need  to  create  a  data  set  file  with  the  same  value  repeated  for  each 
node  in  the  mesh.  If  a  constant  value  is  not  appropriate,  the  ICC  card  specifies 
that  the  initial  condition  varies  spatially  and  the  vdues  are  defined  by  a  data  set 
file. 


Card  Type 

ICC 

Description 

Constant  or  variable  concentration  initial  condition. 

Required 

NO 

Format 

ICC  icon  conval 

Sample 

ICC  0  0.0 

Field 

Variable 

Value 

Description 

1 

ICON 

0 

Constant  concentration  initial  conditions. 

1 

Spatially  variable  concentration  initial  conditions. 

2 

CONVAL 

± 

If  ICON=0,  the  concentration  value.  Othenwise,  omit. 
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Constant  or  variable  head  (ICH) 


The  ICH  is  similar  to  the  ICC  card  in  that  it  can  be  used  as  a  quick  and  easy 
way  to  define  a  cold  start  initial  condition  for  pressure  head.  The  ICH  card  can 
be  used  to  designate  that  the  pressure  head  varies  spatially  and  that  the  values 
will  be  read  from  a  data  set  file.  Alternately,  the  card  can  be  used  to  specify  a 
constant  value  of  total  head.  The  pressure  head  initial  condition  is  then  defined 
in  FEMWATER  by  subtracting  the  elevation  of  each  node  from  the  specified 
total  head  to  define  the  pressure  head  for  the  node. 


Card  Type 

ICH 

Description 

Constant  or  variable  head  initial  condition. 

Repuired 

NO 

Format 

ICC  ihead  hconst 

Sample 

ICC  0100.0 

Field 

Variable 

Value 

Description 

1 

IHEAD 

0 

Constant  total  head  initial  conditions. 

1 

Spatially  variable  pressure  head  initial  conditions. 

2 

HCONST 

± 

If  IHEAD=0,  the  total  head  value.  Otherwise,  omit. 

Initial  condition  file  format  (ICF) 

With  the  exception  of  the  constant  value  approach  for  cold  starts,  the  initial 
conditions  for  cold  starts,  hot  starts,  and  flow  solutions  are  defined  in  data  set 
files.  Data  set  files  can  be  in  either  text  or  binary  format.  The  ICF  card  is  used 
to  designate  the  format  of  the  files  so  that  they  can  be  read  properly  by 
FEMWATER.  A  format  flag  is  specified  for  the  cold  start  or  hot  start  files  and 
another  format  flag  is  specified  for  the  flow  solution  files  used  for  a  transport- 
only  simulation  (if  necessary).  In  addition,  the  card  includes  a  flag  indicating 
whether  the  flow  solution  is  steady  state  or  transient. 


Card  Type 

ICF 

Description 

Initial  condition  file  format. 

Repuired 

YES 

Format 

ICF  icfile  ivfiie  iffu 

Sample 

ICF  1  1  1 

Field 

Variable 

Value 

Description 

1 

ICFILE 

0 

The  initial  condition  files  are  in  text  format. 

1 

The  initial  condition  files  are  in  binary  format. 

2 

IVFILE 

0 

The  flow  files  are  in  text  format. 

1 

The  flow  files  are  in  binary  format. 

3 

IFFU 

0 

Flow  files  contain  a  steady-state  solution. 

1 

Flow  files  contain  a  transient  solution. 
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Super  File  Input 


As  mentioned  in  the  previous  paragraph,  initial  conditions  are  typically 
defined  with  a  set  of  data  set  files.  The  data  set  files  should  be  located  in  the 
same  directory  as  the  other  input  files.  The  names  of  the  data  set  files  to  be  used 
for  initial  conditions  must  be  specified  in  the  FEMWATER  super  file  described 
in  Chapter  2.  The  initial  condition  cards  in  the  super  file  are  shown  in  Figure  9. 
Each  card  consists  of  a  header,  which  identifies  the  type  of  file,  and  the 
filename.  Only  the  cards  corresponding  to  files  needed  for  a  particular 
simulation  should  be  included. 


ICHD  filename 

/*  Pressure  head  init.  cond.  */ 

ICMC  filename 

/*  Moisture  content  (nodal)  init.  cond.  */ 

ICVL  filename 

r  Velocity  init.  cond.*/ 

ICCN  filename 

r  Concentration  init.  cond.  */ 

FLVL  filename 

/*  Velocity  flow  file  (optional)*/ 

FLPH  filename 

/*  Pressure  head  file  (required)*/ 

Figure  9.  The  initial  condition  cards  in  the  super  file 

The  super  file  is  also  used  to  designate  the  names  of  the  solution  files. 
During  a  hot  start,  it  is  possible  to  use  the  same  names  for  the  hot  start  files  and 
the  solution  files.  If  the  same  names  are  used,  FEMWATER  will  append  the 
solution  to  the  end  of  the  hot  start  files  rather  than  overwriting  the  files  or 
creating  new  files. 
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8  Sample  Applications 


The  purpose  of  this  section  is  to  provide  typical  examples  of  FEMWATER 
applications  after  which  one  might  pattern  new  applications.  Five  examples  are 
provided  that  demonstrate  different  capabilities  of  FEMWATER.  Each 
discusses  boundary  and  initial  condition  specification.  The  input  files  associated 
with  these  examples  are  distributed  with  the  FEMWATER  executable  and  can  be 
found  in  the  “examples”  directory. 


Problem  1 :  Steady-State  Wellhead  Protection 

The  initial  support  for  developing  FEMWATER  within  the  GMS  was 
provided  by  the  AERL  to  advance  the  state  of  the  art  and  provide  practice  in 
doing  wellhead  protection  studies.  As  a  result,  of  the  technologies  accepted  by 
the  U.S.  Environmental  Protection  Agency  (EPA)  to  do  wellhead  protection 
studies,  FEMWATER  and  GMS  provide  Ae  most  rigorous  modeling  tools. 
Several  of  these  studies  have  been  completed  within  the  Department  of  Defense 
and  many  more  are  in  the  planning  phase. 

The  provided  example  illustrates  how  to  use  FEMWATER  to  conduct  a 
steady-flow  simulation  for  a  wellhead  protection  problem.  The  left  and  right 
sides  of  the  model  are  bounded  by  freshwater  head  boundary  conditions  (Figure 
10).  The  front  and  back  of  the  model  are  groundwater  divides  represented  by  no¬ 
flux  boundary  conditions.  The  bottom  of  the  model  is  an  impervious  stratum 
that  does  not  allow  significant  amounts  of  water  to  leave  the  model.  This  is 
commonly  located  at  aquicludes  or  deep  aquitards.  The  pressure  head  is 
assumed  to  be  hydrostatic  throughout  the  model  with  water  table  initial 
conditions  assumed  to  be  a  constant  12.0  m.  The  following  boundary  conditions 
are  given:  constant  head  of  12.0  m  on  the  left  side  of  the  model,  13.0  m  on  the 
right  side  of  the  model,  and  no  flux  on  the  front  and  back  of  the  model. 

Seventeen  pumping  wells  are  located  in  the  top  layer  of  coarse  sand  near  the 
middle  of  the  region.  The  pumping  rate  of  the  wells  is  a  constant  0.25  m^/hr. 

One  of  the  wells  could  be  increased  to  determine  the  effect  on  the  other  wells. 
The  head  tolerance  is  0.001  m.  The  pointwise  iterative  solver  is  used  to  solve 
the  matrix  equations.  The  computational  mesh  consists  of  4,000  elements  and 
3,507  nodes. 
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Figure  10.  Wellhead  protection  application  showing  numerical  model  mesh  and 
assigned  boundary  conditions 

This  problem  took  approximately  1.1  minutes  for  the  steady-state  solution  to 
converge  on  a  200-MHz  Pentium  class  personal  computer  with  32  MB  of 
memory. 


Problem  2:  Transient  Confined  Disposal 

Confined  Disposal  Facilities  (CDF’s)  are  essentially  settling  basins  that 
accept  dredged  materials  from  navigation  chaimels.  The  dredged  material 
typically  enters  the  CDF’s  as  aqueous  slurries  for  the  purposes  of  settling  out  the 
solids  and  returning  the  clean  water  to  the  sea.  This  is  a  common  activity  for  the 
U.S.  Army  Corps  of  Engineers.  Regulatory  agencies  are  increasingly  requiring 
the  use  of  3-D  groundwater  models  to  determine  the  path  of  groundwater  from 
the  CDF  to  surrounding  areas.  Since  dredged  material  sometimes  contains 
compounds  that  should  not  enter  drinking  water  supplies,  models  of  CDF’s  can 
determine  if  the  contamination  of  drinking  water  wells  was  caused  by  a 
particular  CDF.  In  the  interest  of  prevention,  CDF’s  can  be  sited  and  designed  to 
ensure  that  they  will  not  contaminate  local  wells. 

The  provided  example  is  presented  to  illustrate  the  procedure  for  conducting 
a  transient  flow  simulation  at  a  CDF.  The  left  and  right  sides  of  the  model  are 
bounded  by  freshwater  head  boundary  conditions  (Figure  11).  The  front  and 
back  of  the  model  are  bounded  by  groundwater  divides  represented  by  no-flux 
boundary  conditions.  The  bottom  of  the  model  is  an  impervious  stratum.  A 
pond  (CDF)  filled  with  a  conservative  chemical  constituent  is  located  on  the 
right  side  of  the  model.  Pressure  is  assumed  to  be  hydrostatic  throughout  the 
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model  with  water  table  initial  conditions  assumed  to  be  a  constant  12.0  m.  The 
following  boundary  conditions  are  given:  constant  head  of  12.0  m  on  the  left  side 
of  the  model,  variable  head  with  slowly  rising  water  table  assumed  on  the  right 
side  of  the  model,  constant  head  of  13.0  m  with  concentration  of  10  mg/^ 
assumed  in  the  pond,  and  no  flux  imposed  on  the  front  and  back  of  the  model.  A 
well  with  multiple  screen  openings  is  located  near  the  center  of  the  region.  The 
pumping  rate  of  the  well  is  varied  slowly  with  a  slightly  increasing  rate.  It  starts 
with  0.05  m^/hr  and  increased  to  0.1  m^/hr  at  time  =  10.0  hr,  where  it  stays 
constant  at  0.1  m%r.  The  head  tolerance  is  0.001  m  and  the  pointwise  iterative 
solver  is  used.  A  total  of  200  hr  with  variable  time-step  is  simulated. 


Figure  11.  Confined  disposal  facility  application  showing  numerical  model 
mesh  and  assigned  boundary  conditions 

This  problem  took  approximately  9.1  minutes  for  the  transient-flow  solution 
to  converge  on  a  200-MHz  Pentium  class  personal  computer  with  32  MB  of 
memory. 


Problem  3:  Transient  Groundwater  Remediation 

This  example  illustrates  a  means  for  modeling  a  simple  remediation  situation 
using  a  cutoff  wall  to  capture  flows  from  a  chemical-laden  pond.  The  left  and 
right  sides  of  the  model  are  bounded  by  freshwater  head  boundary  conditions 
(Figure  12).  The  front  and  back  of  the  model  are  bounded  by  groundwater 
divides  represented  by  no-flux  boundary  conditions.  The  bottom  of  the  model  is 
an  impervious  stratum.  A  cutoff  wall  consisting  of  impervious  media  is  installed 
in  the  middle  of  the  model  and  a  chemical  lagoon  is  located  on  the  right  side. 
Pressure  head  is  assumed  hydrostatic  throughout  the  region  with  water  table 
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initial  conditions  of  12.0  m  throughout.  The  boundary  conditions  are  as  follows: 
constant  head  of  12.0  m  on  the  left  side,  variable  head  with  slowly  rising  head  on 
the  right  side,  constant  head  of  13.0  m  with  a  concentration  of  10  mg/ i  in  the 
ponding  lagoon,  and  no-flux  conditions  on  the  front,  back,  and  bottom.  The 
pumping  wells  are  located  in  the  coarse  sand  layer  near  the  center  of  the  region. 
The  pumping  rate  of  the  well  slowly  increases.  The  pumping  starts  with  0.05 
m^/hr  and  increases  to  0.1  m^/hr  at  time  =  10.0  hr,  where  it  stays  constant  at  0.1 
m^/hr.  The  pointwise  iterative  solver  is  used.  A  total  of  500  hr  with  variable 
time-step  is  simulated. 


Figure  12.  Groundwater  remediation  application  showing  numerical  model  mesh 
and  assigned  boundary  conditions 

This  simulation  took  approximately  14.2  minutes  for  the  flow  and  transport 
solution  to  converge  on  a  200-MHz  Pentium  class  personal  computer  with  32 
MB  of  memory. 


Problem  4:  Transient  Chemical  Spill 

This  example  presents  an  example  of  how  to  model  a  chemical  spill  under 
transient  conditions.  The  left  and  right  sides  are  bounded  by  freshwater  constant 
head  boundary  conditions  (Figure  13).  The  front  and  back  of  the  region  are 
bounded  by  groundwater  divides  represented  by  no-flux  boundary  conditions. 
The  bottom  of  the  model  is  bounded  by  an  impervious  stratum.  The  chemical 
spill  is  located  near  the  center  of  the  region.  The  pressure  is  assumed  to  be 
hydrostatic.  The  initial  water  surface  elevation  is  a  constant  1 1.0  m.  The 
following  boundary  conditions  are  given:  a  constant  head  of  1 1.0  m  on  the  left 
side  of  the  region,  a  variable  head  with  a  slowly  rising  water  table  on  the  right 
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side  of  the  region,  a  liquid  chemical  with  a  concentration  of  10.0  mg/  i  spilled 
on  the  ground  surface,  rainfall  imposed  as  a  variable  boundary  condition  on  the 
ground  surface,  and  no  flux  imposed  on  the  front  and  back  of  the  region.  A  set 
of  wells  are  located  near  the  center  of  the  region.  The  pumping  rate  of  the  well 
is  varied  with  a  slowly  increasing  rate.  It  starts  with  0.05  m^/hr  and  increases  to 
0.1  m^/hr  at  time  =  10.0  hr,  and  it  stays  constant  at  0.1  m^r.  The  head  tolerance 
is  0.001  m  and  the  pointwise  iterative  solver  is  used.  A  total  of  150  hr  with 
variable  time-steps  is  simulated. 


Figure  13.  Chemical  spill  application  showing  numerical  model  mesh  and 
assigned  boundary  conditions 

This  problem  took  approximately  17.8  minutes  for  the  flow  and  transport 
solution  to  converge  on  a  200-MHz  Pentium  class  personal  computer  with  32 
MB  of  memory. 


Problem  5:  Transient  Salinity  Intrusion 

The  salinity  intrusion  example  demonstrates  how  FEMWATER  can  be  used 
to  study  a  variety  of  coastal  aquifer  problems.  It  is  currently  being  used  by  the 
Corps  of  Engineers  to  determine  the  impacts  of  deepened  navigation  channels  on 
salinity  intrusion  in  coastal  aquifers  and  how  increased  groundwater  pumping 
will  affect  drinking  water  supplies. 

In  order  to  correctly  simulate  salinity  movement  in  coastal  aquifers,  it  is 
necessary  to  have  coupled  flow  and  density-driven  transport.  The  example 
problem  is  a  schematic  example  of  a  coastal  salinity  intmsion  problem.  The  left 
side  of  the  region  is  bounded  by  seawater,  and  the  right  side  is  bounded  by  fresh 
water  (Figure  14).  The  front  and  back  of  the  region  are  bounded  by  groundwater 
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divides  and  the  bottom  by  an  impervious  stratum.  A  pumping  well  is  located 
near  the  center  of  the  region.  The  boundary  conditions  are  specified  as  a  20.0- 
kg/^  tidal  boundary  on  the  left  side  and  fresh  water  with  a  rising  water  table  on 
the  right  side.  The  pumping  rate  of  the  well  is  varied  with  a  slowly  increasing 
rate.  It  starts  with  0.5  nr/hr  and  increases  to  2.0  m^/hr  at  time  =  100.0  hr  where 
it  stays  constant  at  2.0  m^/hr.  A  total  of  10.0  hr  with  variable  time-steps  is 
simulated.  The  computational  mesh  consists  of  5,664  elements  and  3,484  nodes. 


Figure  14.  Transient  salinity  Intmsion  application  showing  numerical  model 
mesh  and  assigned  boundary  conditions 

This  problem  took  approximately  45.7  minutes  for  the  transient  solution  to 
converge  on  a  200-MHz  Pentium  class  personal  computer  with  32  MB  of 
memory. 
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Appendix  A  Mathematical 
Formulation 


Governing  Equations  for  Flow 


From  the  notion  for  continuity  of  fluid,  continuity  of  solid,  consolidation  of 
the  media,  and  the  equation  of  state  (Yeh  1992)'  are  obtained  the  starting 
equation  for  this  derivation: 


■(Vp  +  pgVz) 


-V.(pn,SV,)+p*q  = 


^(neSp) 

at 


(Al) 


where 

p  =  fluid  density  (M/L^) 

k  =  intrinsic  permeability  tensor  of  the  media  (L^) 

p  =  dynamic  viscosity  of  the  fluid  (M/L/T) 

p  =  fluid  pressure  [(ML/T^)/L^] 

g  =  acceleration  of  gravity  (L/T^) 

z  =  potential  head  (L) 

rie  =  effective  porosity  (L^/L^) 

S  =  degree  of  samration  (dimensionless) 

Vj  =  velocity  of  the  deformable  surface  due  to  consolidation  (L/T) 


‘  References  cited  in  this  Appendix  are  located  at  the  end  of  the  main  text. 
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p*  =  density  of  the  injected  fluid  (M/L^) 
q  =  internal  source/sink  [(L^/T)/LL^] 
t  =  time  (T) 

Expanding  the  right  hand  of  Equation  A1 : 


^(PeSp) 

at 


as 


Expanding  Equation  A2  by  the  chain  rule: 


a(neSp) 

at 


_c„ 

^^apat 


+  Sng 


ap  ac 
ac  at 


+ps- 


an. 


as 

■z-^  +  n.p— 

at  "'^at 


(A2) 


(A3) 


where  C  is  chemical  concentration  (M/L^).  Rearranging  Equation  A3  gives: 


^(^eSp)  ^iP, 

at  ^  ^apat  P^at 


ap  ac 


as 

at 


(A4) 


where  the  first  and  second  terms  represent  the  storativity  term,  the  third  term  is 
the  density-concentration  coupling  term,  and  the  fourth  term  is  the  unsaturated 
term.  Substituting  Equation  A4  into  Equation  Al: 


V- 


^•(Vp  +  pgVz) 

M' 


+p  q 


ap  ap  ap  ac  as 

d^Ji  dc^ 

3  n 

+  ps-^ +SpV .  n,  V.  +  n.  V.  •  V{Sp) 


(A5) 


Making  the  approximation  by  neglecting  the  second-order  term: 


n,V,-V(Sp)-- - >0 


(A6) 


gives 
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V- 


^.(Vp  +  pgVz) 


+p  q  = 


3p  0p  3p  9C  dS 
+  pS^+SpV-n,V, 


(A7) 


Defining  the  compressibility  of  the  fluid  as: 


(A8) 


where  P  is  the  compressibility  of  the  fluid  (LT^/M).  Also  defining  the  moisture 
content  as: 


e  =  sn. 


(A9) 


where  6  is  the  moisture  content  (dimensionless).  Substitute  Equations  A8  and 
A9  into  Equation  A7  and  rewrite  it  to  obtain: 


V 


^■(Vp+pgVz) 


^  3p  9p  8C 

epp^+e 


+  p  q  = 
9S 


9c  at  "^“^^at 


+  pS 


at 


V(n,V.) 


(AlO) 


Remembering  the  continuity  statement  of  incompressible  solids,  a 
compressible  skeleton  is  defined  as  (Yeh  1992): 


a(i-nj 

at 


+  V-(l-nJV,=0 


(All) 


Rearranging  Equation  A1 1  in  the  following  form: 

^+v.„.v.=v.v. 


(A12) 
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Substituting  Equation  A12  into  Equation  AlO  gives: 


V- 


^^■{Vp  +  pgVz) 
L 


+  p  q  = 


dp  d p  dC  d S 


(A13) 


Recalling  that  the  flux  of  solid  velocity  is  the  divergence  of  Fj  (Yeh  1992): 
V-V,=a-^  (A14) 


where  a  is  the  coefficient  of  consolidation  of  the  media  (LT^/M).  Substituting 
Equation  A14  into  Equation  A13  and  rewriting: 


V- 


^■(Vp  +  pgVz) 


+  p  q 


9p  3p  dC  9S 
P(ep  +  Sa) —  +  0^-^  +  Hep-^ 


(A15) 


Substituting  Equation  A9  into  Equation  15  gives 


V- 
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^•(Vp  +  pgVz) 


+  p  q 


0 


0P  + — a 


V 


n 


e  y 


9p  9p  dC  ds 

9t  ac  9t  "^at 


(A16) 


Experimental  evidence  has  shown  that  the  degree  of  saturation  is  a  function  of 
pressure  as: 

S  =  S(p)  (A17) 

Substitution  of  Equation  A17  into  Equation  A16  gives: 
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'•  ^(Vp  +  pgVz) 


+  p  q  = 


^  0  9p  9C  dS  9p 

p  0B  + — a  —  +  0-,^-,  +pne. 

M  ^  n,  J9t  9C  at  dp  at 


(A18) 


Next,  the  reference  pressure  head  is  defined  as: 


(A19) 


where  h  is  the  reference  pressure  head  (L)  and  po  is  the  reference  water  density 
(M/L^). 


Substituting  Equation  A19  into  Equation  A18  gives: 

V-  ^-{Pog^h  +  pgVz)  +p*q  = 

^  0  ^  dh  dpdC  dS  dh 

p  68  + — a  +  +pn.  „ 

I  n.  r  dt  dC  dt  ‘  dh  dt 


(A20) 


Dividing  Equation  A20  by  Po  and  rearranging  yields: 

V.  Vzl  +—q  = 

_  P  \  Po  Po 

p(  6  ']dh  6  dpdC  p  dSdh 

— I  1 57 + ^^57 


(All) 


Defining  the  modified  compressibilities  of  the  media  and  water  as 


a  =apog 


(A22) 


P'  =  PPo8 


(A23) 


where  a’  is  the  modified  compressibility  of  the  media  (1/L)  and  P'  is  the 
modified  compressibility  of  the  water  (1/L).  Substituting  Equations  A22  and 
A23  into  Equation  A21  and  rearranging: 
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V-  — •[vh  +— Vz]  +— q  = 

I  \^  {  Po  jj  Po 

p  f  ,  6  ds"iah  e  ap  ac 

oc  +B  0  +  n.“~  "T  + 

Pol  n,  dhJat  po  ac  at 


(A24) 


Defining  the  storage  coefficient  as: 


e  dS 

F  =  a' —  +  P'0  +  ne  — 
n„  ®  dh 


(A25) 


where  F  is  the  storage  coefficient.  Substituting  Equation  A25  into  Equation  A24 
and  following  Frind  (1982)  by  neglecting  the  second  term  on  the  right  side  of 
Equation  A24, 


V.  M.fvh+^Vz 

R  I  Po  j  Po  Po 


(A26) 


Defining  the  relation: 


(A27) 


where  iT  is  the  hydraulic  conductivity  tensor.  Substituting  Equation  A27  into 
Equation  A26  and  rearranging  gives  the  density-dependent  flow  equation: 


— f|^  =  V-  K-fvh-H-^Vz]  -h— q 

Po  dt  I  Po  j  Po 


(A28) 


From  Darcy’s  law: 

V  =  ■  (Vp  +  pgVz) 

p  |X  '  ' 


(A29) 


where  V  is  the  Darcy  flux  (L/T).  Recalling  Equation  A19  and  substituting  into 
Equation  A29  gives: 


V  =  •  (pogVh  -b  pgVz) 

P  M- 


(A30) 
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Rearranging  Equation  A30: 


V  =  - 


pgk 


^Vh  +  Vz 

Ip  ; 


(A31) 


and  substituting  Equation  A27  into  Equation  A31,  we  get  the  Darcy  flux 
equation  for  density-dependent  flow  in  its  final  form: 


V  = 


-K 


^Vh-f-Vz 

VP 


(A32) 


The  density  and  dynamic  viscosity  are  functions  of  chemical  concentration 
and  are  assumed  to  take  the  following  form: 


-^  =  a, -HajC-i-ajC'+a^C^  (A33) 

Po 


and 


=  aj  +  agC  +  SiyC^  +  agC^ 

Po 


(A34) 


where  C  is  the  chemical  concentration  (M/L^)  and  aj,  02, aj,  as  are  the 
parameters  (L^/M)  that  are  used  to  describe  the  concentration  dependence  of 
water  density  and  dynamic  viscosity. 

In  the  specific  case  of  seawater  intrusion,  the  constitutive  relation  between 
fluid  density  and  concentration  takes  the  following  form: 

p  =  p„(l-l-ec)  (A35) 


where  c  is  the  dimensionless  chemical  concentration  (the  actual  one  divided  by 
the  maximum  one)  and  e  is  the  dimensionless  density  reference  ratio  defined  as: 


e  = 


(A36) 


where  pmax  is  the  maximum  density  of  the  fluid  (M/L^)  and  Po  is  the  reference 
(freshwater)  density  (M/L^).  It  should  be  noted  that  Equation  A36  implicitly 
assumes  that  the  fluid  is  incompressible  and  under  isothermal  conditions 
(Galeati,  Gambolati,  and  Neumann,  1992). 
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The  initial  conditions  for  the  flow  equations  are  stated  as: 


h  =  h;  (x,  z)  in  R 


(A37) 


where  R  is  the  region  of  interest  and  hi  is  the  prescribed  initial  condition  for 
hydraulic  head.  The  specification  of  boundary  conditions  is  probably  the  most 
critical  and  complex  chore  in  flow  modeling.  As  explained  by  Yeh  (1987),  the 
boundary  conditions  of  the  region  of  interest  can  be  examined  from  a  dynamic, 
physical,  or  mathematical  point  of  view.  From  a  dynamic  standpoint,  a  boundary 
segment  can  be  considered  either  as  impermeable  or  flow-through.  On  the  other 
hand,  from  a  physical  point  of  view,  such  a  segment  could  be  classified  as  a  soil- 
soil  interface,  soil-air  interface,  or  soil-water  interface.  Lastly,  from  a 
mathematical  point  of  view,  the  boundary  segment  can  be  classified  as  one  of 
four  types  of  boundary  conditions,  namely  as  (a)  Dirichlet,  (b)  Neumann,  (c) 
Cauchy,  or  (d)  variable  boundary  conditions.  In  addition,  a  good  numerical 
model  must  be  able  to  handle  these  boundary  conditions  when  they  vary  on  the 
boundary  and  are  either  abruptly  or  gradually  time-dependent. 

The  Dirichlet  boundary  condition  is  usually  applied  to  soil-water  interfaces, 
such  as  streams,  artificial  impoundments,  and  coastal  lines,  and  involves 
prescribing  the  functional  value  on  the  boundary.  The  Neumann  boundary 
condition,  on  the  other  hand,  involves  prescribing  the  gradient  of  the  function  on 
the  boundary  and  does  not  occur  very  often  in  real-world  problems.  This 
condition,  however,  can  be  encountered  at  the  base  of  the  media  where  natural 
drainage  occurs.  The  third  type  of  boundary  condition,  the  Cauchy  boundary 
condition,  involves  prescribing  the  total  normal  flux  due  to  the  gradient  on  the 
boundary.  Usually  surface  water  bodies  with  known  infiltration  rates  through 
the  bottom  layers  of  their  sediments  or  liners  into  the  subsurface  media  are 
administered  this  boundary  condition.  If  a  sod-air  interface  exists  in  the  region 
of  interest,  a  variable  boundary  condition  is  employed.  In  such  a  case,  either 
Dirichlet  or  Cauchy  boundary  conditions  dominate,  mainly  depending  on  the 
potential  evaporation,  the  conductivity  of  the  media,  and  die  availability  of  water 
such  as  rainfall  (Yeh  1987a). 

From  this  discussion,  four  types  of  boundary  conditions  can  be  specified  for 
the  flow  equations  depending  on  the  physical  location  of  the  boundaries: 


a.  Dirichlet  boundary  conditions: 

h  =  h^{x^,y^,Zh,t)  onBj,  (A38) 

b.  Neumarm  boundary  conditions: 
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-n  K 


^■Vh 

VP  ) 


(A39) 


c.  Cauchy  boundary  conditions: 


n  K 


Po 


V/i  + Vz 


■■qcix^^yi,,Zi,t)  onB^, 


(A40) 


d.  Variable  boundary  conditions  -  during  precipitation  period: 


h  =  hp(xi,,y^,Zh,t)  on  B^, 


(A41) 


or 


-n  K 


fp  ^ 

^Vh  +  Vz 

VP  J 


=  qp{xb,yb>Zb,t)  on  B^,  (A42) 


e.  Variable  boundary  conditions  -  during  nonprecipitation  period: 


h  =  hp(xb,yb,Zb,t)  on  B,, 


(A43) 


or 


h  =  h„,(xb,yb,Zb,t)  on  B,, 


(A44) 


or 


-n  K 


— Vh  +  Vz 

VP  y 


=  qe(xb.yb’Zb’t)  on  B,,  (A45) 


f.  River  boundary  conditions: 
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-n  K 


^Vh  +  Vz 

.P  J 


=  ^(hR-h)  on  B,, 


(A46) 


where 

n  =  outward  unit  vector  normal  to  the  boundary 

(xb,  yh  Zb)  =  spatial  coordinate  on  the  boundary 

hd  =  Dirichlet  functional  value 

qh  =  Neumann  flux 

qc  =  Cauchy  flux 

=  Dirichlet  boundary 

Bn  =  Neumann  boundary 

Be  =  Cauchy  boundary 

fiv  =  variable  boundary 

hp  =  ponding  depth  on  the  variable  boundary 

qp  =  throughfall  of  precipitation  on  the  variable  boundary 

hm  =  minimum  pressure  on  the  variable  boundary 

qe  =  evaporation  rate  (potential  evaporation)  on  the  variable  boundary 

Kr  =  hydraulic  conductivity  of  the  river  bottom  sediment  layer 

bR  =  thickness  of  the  river  bottom  sediment  layer 

Hr  =  depth  of  the  river  bottom  measured  from  the  river  surface. 

Note  that  only  one  of  Equations  A41-A45  is  utilized  at  any  point  on  the  variable 
boundary  at  any  time. 


Governing  Equations  for  Transport 

The  governing  equations  for  material  transport  through  groundwater  systems 
are  derived  based  on  the  laws  of  continuity  of  mass  and  flux.  The  major 
processes  that  are  included  are  advection,  dispersion/diffusion,  decay, 
adsorption,  biodegradation  through  both  liquid  and  solid  phases,  the 
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compressibility  of  media,  as  well  as  source(s)/siiik(s).  Let  C  be  the  dissolved 
concentration  and  S  be  the  adsorbed  concentration.  The  governing  equation  of 
the  spatial-temporal  distribution  of  dissolved  concentrations  can  be  obtained  by 
applying  this  principle  of  mass  balance  in  integral  form  as  follows: 

^  J  (0  C  p^s)dv  =  -\n{e  c)Vfe  t/r 

V  r 

-jn-Jdr-l{eK^C  +  p,K^s)dv  (A47) 

r  V 

- 1  X{dC  -1-  Pi,S)dv  -l-| mdv 

V  V 

where 

V  =  material  volume  containing  constant  amount  of  media  (L^) 

C  =  dissolved  concentration  (M/L^) 

S  =  adsorbed  concentration  (M/M) 
ph  =  bulk  density  of  the  medium  (M/L^) 
r=  surface  enclosing  the  material  volume  v  (L^) 
n  =  outward  unit  vector  normal  to  the  surface  G 
Vfs  =  fluid  velocity  relative  to  the  solid  (L/T) 

J  =  surface  flux  with  respect  to  fluid  velocity  [(M/T)/L^] 

Kw  =  first-order  biodegradation  rate  constant  through  dissolved  phase  (1/T) 

Ks  =  first-order  biodegradation  rate  constant  through  adsorbed  phase  (1/T) 

X  =  decay  constant  (1/T) 

m  =  external  source/sink  rate  per  medium  volume  [(M/L^)/T] 

By  the  Reynolds  transport  theorem  (Owczarek  1964),  Equation  A47  can  be 
written  as 

J  dv + J  n .  (e  c + P  ,s)  V.  dr + J  n .  (e  cv„ )  dr = 

V  T  T  (A48) 

-  jn  ■  J  dr  -J  (bK^C  -t-  PbK3S)dv  -  j  X(e  C  +  p,s)dv  -t-  J  mdv 

r  V  V  V 
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where  is  the  velocity  of  the  solid.  The  fluid  velocity  relative  to  the  solid 
and  Darcy  velocity  V  are  related  to  each  other  by 


V  =  0%, 


(A49) 


Applying  the  Gaussian  divergence  theorem  to  Equation  A48  and  using  the 
fact  that  V  is  arbitrary,  one  can  obtain  the  following  continuity  in  differential 
form: 

9(ec+p,s)^^  [v, (ec + p,5)] + V ■  (vc) = -v ■  j 

-  {SK^C + p,  X,s)  -  A(«C + p,5)  +  m 


The  second  term  of  Equation  A50  can  be  expressed  as 

V .  [v.  (ec + P,s)] = V(9C + P,s)  ■  V, + (ec + P,s)v  •  V.  (as  d 


The  first  term  on  the  right  side  of  Equation  A5 1  is  the  product  of  two  small 
vectors  and  will  be  neglected.  If  aU  the  displacement  of  the  medium  is  assumed 
to  be  vertical  (e.g.,  vertical  consolidation),  the  solid  velocity  becomes 


V-V, 


(A52) 


The  surface  flux  J  has  been  postulated  to  be  proportional  to  the  gradient  of  C 
(Nguyen  et  al.  1982)  as 

J  =  -0D-VC  (A53) 

GD  =  a^ I  V|6  +  (ul  -  )  W  / 1  V|  +  a^et5  (A54) 

where 

ot  =  transverse  diffusivity  (L) 

8  =  Kronecker  delta  tensor 
I  FI  =  magnitude  of  the  Darcy  velocity  V  (L/T) 
ai  =  longitudinal  diffusivity  (L) 
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am  =  molecular  diffusion  coefficient  (L^/T) 

T=  tortuosity 

Substitution  of  Equations  A51  through  A54  into  Equation  A50  yields 


3(ec+p.s) 


3t 


(X":: — 1-  X 

3t  j 


+v-(vc)-v-(eD-vc)  = 

|(9C  +  p.s)-(eK.C  +  p.K.s)  +  m 


(A55) 


Equation  A55  is  simply  the  statement  of  mass  balance  over  a  differential  volume. 
The  first  term  represents  the  rate  of  mass  accumulation,  the  second  term 
represents  the  net  rate  of  mass  flux  due  to  advection,  the  third  term  is  the  net 
mass  flux  due  to  dispersion  and  diffusion,  the  fourth  term  is  the  rate  of  mass 
production  and  reduction  due  to  consolidation  and  radioactive  decay,  the  fifth 
term  is  the  degradation  rate  due  to  microbial  transformation  through  aqueous  and 
adsorbed  phases,  and  the  last  term  is  a  source/sink  term  corresponding  to 
artificial  injection  and  or  withdrawal. 


Equation  A55  is  written  in  conservative  form.  Using  the  advective  form  is 
sometimes  more  appropriate,  especially  if  the  finite  element  method  is  used  to 
simulate  the  chemical  transport  equation.  More  importantly,  an  advective  form 
of  the  transport  equations  allows  one  to  use  the  mixed  Lagrangian-Eulerian 
approach,  which  can  better  solve  advection-dominant  transport  problems  (Yeh 
and  Tripathi  1987).  Therefore,  an  advective  form  of  the  transport  equation  is 
derived  by  expanding  the  advection  term  and  using  the  continuity  equation  for 
water  flow.  The  water  flow  equation  can  be  rewritten  in  the  following  form: 


p  9h 
Po 


=  -V- 


^po  ; 


v-v-v-v 

Po  Vpo / 


(A56) 


which  is  conservation  of  fluid  mass.  Substituting  Equation  A56  into  A55  and 
performing  the  necessary  manipulation  gives 
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e^+Pb|^+v-vc-v-(0DVc)  = 

3t  dt 


'  Kec+p»s)-{9K.c+p.K.s)+ 


a'-z— +X 


V 


at 


(A57) 


m- 


» 

p 

ah  p. ^ 

fp^ 

ae' 

— qC  + 
P 

F— +— V-V 

1  at  p 

vPo  J 

"atj 

Equation  A57  involves  two  unknowns  C  and  5,  so  constitutive  relationships 
must  be  posed.  For  the  present  model,  the  following  empirical  relationships  are 
used: 


S  =  K,C 

for  linear  isotherm 

(A58) 

Cl 

for  Langmuir  isotherm 

(A59) 

“  l  +  KC 

S  =  KC" 

for  Freundlich  isotherm 

(A60) 

where 

Ka  =  distribution  coefficient  (L^/M) 

Smax  =  maximum  concentration  of  the  medium  in  the  Langmuir  nonlinear 
isotherm 

K  =  coefficient  in  the  Langmuir  or  Freundlich  nonlinear  isotherm 

n  =  power  index  in  the  Freundlich  nonlinear  isotherm. 

In  order  to  use  the  Lagrangian-Eulerian  approach,  Equation  A57  is  rewritten 
as: 

a.  For  Linear  Isotherm  model: 
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(0  +  A*^.)^  =  V(eDVC)- 

* 

(eKS  +  P,K,S)  +  m-^qC^ 


( 


,dh 


A 


OC'-r — ^•  A 

V  y 


V 

(p"' 

r  "T — 1 - V  *  V 

[  dt  p 

<PoJ 

dt  ^ 

(ec+Pi,^)- 


(A61) 


v;  = 


0+PbKd 


(A62) 


b.  For  Langmuir  and  Freundlich  models : 


C  dS  8C  /  X 


cx  "T — + A. 

at 


l(ec+p,s)-(eK,c+p,K.s)+ 


m- — qC  + 
P 


dh  P-  __  _ 

(P) 

F— +— V-V 

t  dt  p 

<Po, 

ae 

at 


where 


dJ) 

Dt 


and 


D.,() 

Dt 


(A63) 


(A64) 


denote  the  material  derivatives  of  ( )  with  respect  to  time  in  reference  to  the 
retarded  velocity  Vd  and  fluid  velocity  Vf,  respectively. 

To  complete  the  description  of  the  hydrological  transport  as  given  by 
Equations  A61  or  A63,  initial  and  boundary  conditions  must  be  specified  in 
accordance  with  dynamic  and  physical  considerations.  It  will  be  assumed  that 
initially  the  dissolved  concentrations  are  known  throughout  the  region  of 
interest,  that  is, 
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C  =  Cj  (x,  z)  in  R 


(A65) 


where  C,  is  the  initial  concentration  and  R  is  the  region  of  interest.  Initial 
concentrations  for  the  dissolved  concentrations  may  be  obtained  from  field 
measurements  or  by  solving  the  steady-state  version  of  Equation  A61  or  A63 
with  time-invariant  boundary  conditions. 

The  specification  of  boundary  conditions  is  a  difficult  and  intricate  task  in 
transport  modeling.  From  the  dynamic  point  of  view,  a  boundary  segment  may 
be  classified  as  either  flow-through  or  impervious.  From  a  physical  point  of 
view,  it  is  a  soil-air  interface,  soil-soil  interface,  or  soil-water  interface.  From 
the  mathematical  point  of  view,  it  may  be  treated  as  a  Dirichlet  boundary  on 
which  the  total  analytical  concentration  is  prescribed,  a  Neumann  boundary  on 
which  the  flux  due  to  the  gradient  of  total  analytical  concentration  is  known,  or  a 
Cauchy  boundary  on  which  the  total  flux  is  given.  An  even  more  difficult 
mathematical  boundary  is  the  variable  boundary  conditions  on  which  the 
boundary  conditions  are  not  known  a  priori  but  are  themselves  the  solution  to  be 
sought.  In  other  words,  on  the  mathematically  variable  boundary,  either 
Neumann  or  Cauchy  conditions  may  prevail  and  change  with  time.  Which 
condition  prevails  at  a  particular  time  can  be  determined  only  in  the  cyclic 
processes  of  solving  the  governing  equations  (Freeze  1972a,  1972b;  Yeh  and 
Ward  1980;  Yeh  and  Ward  1981). 

Whatever  point  of  view  is  chosen,  all  boundary  conditions  eventually  must 
be  transformed  into  mathematical  equations  for  quantitative  simulations.  Thus, 
we  will  specify  the  boundary  conditions  from  the  mathematical  point  of  view  in 
concert  with  dynamic  and  physical  considerations.  The  boundary  conditions 
imposed  on  any  segment  of  the  boundary  are  taken  to  be  either  Dirichlet, 
Neumann,  Cauchy,  or  variable.  Thus,  the  global  boundary  may  be  split  into  four 
parts,  Bd,  Bn,  Be,  and  B„,  denoting  Dirichlet,  Neumann,  Cauchy,  and  variable 
boundaries,  respectively.  The  conditions  imposed  on  the  first  three  types  of 
boundaries  are  given  as: 

a.  Prescribed  concentration  (Dirichlet)  boundary  conditions: 

C  =  Cd(Xb,yb,Zb,t)  on  B<j  (A66) 


b.  Neumann  boundary  conditions: 

n-(-eD-Vc)  =  qb(Xb,yb,Zb,t)  on  B„  (A67) 


c.  Cauchy  boundary  conditions: 
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n-(vC-eD-VC)  =  q,(Xb,yb,Zb,t)  on  B 


(A68) 


where 


Cd  =  concentration  on  the  Dirichlet  boundary  Bd 

{Xb,yb,Zb)  =  spatial  coordinate  on  the  boundary 

n  =  outward  unit  vector  normal  to  the  boundary 

q„  =  gradient  flux  through  the  Neumann  boundary 

qc  =  total  flux  through  the  Cauchy  boundary  Be 

The  conditions  imposed  on  the  variable-type  boundary,  which  is  normally 
the  soil-air  interface  or  soil-water  interface,  are  either  the  Neumaim  with  zero 
gradient  flux  or  the  Cauchy  with  given  total  flux.  The  former  is  specified  when 
the  water  flow  is  directed  out  of  the  region  from  the  faraway  boundary,  whereas 
the  latter  is  specified  when  the  water  flow  is  directed  into  the  region.  This  type 
of  variable  condition  would  normally  occur  at  flow-through  boundaries.  Written 
mathematically,  the  variable  boundary  condition  is  given  by 

Ii-(VC-0D-Vc)  =  ii- VCy(Xb,yb,Zi,,t)  onB^  if  n- V  <0 

;  ,  (A69) 

n(-0DVC)  =  OonB,  ifn-V>0 

where  Cy  is  the  specified  concentration  of  water  through  the  variable  boundary 
and  By  is  the  variable  boundary. 
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Appendix  B  Numerical 
Formulation 


Introduction 

The  initial-boundary  value  problem  described  by  the  governing  equations  of 
the  flow  and  transport  modules  of  FEMWATER  along  with  the  boundary 
conditions  cannot,  in  general,  be  solved  analytically  using  applied  mathematics. 
Hence,  in  order  to  solve  these  sets  of  governing  equations,  numerical  methods 
are  the  only  mathematical  tools  capable  of  handling  this  task.  Although  there  are 
many  different  numerical  approximation  methods  capable  of  reducing  partial 
differential  equations  to  simpler  systems  of  algebraic  equations,  there  are  only 
two  numerical  methods  that  are  most  common  and  that  can  be  employed  to  solve 
the  most  basic  form  of  the  governing  equations.  These  are  the  finite  difference 
and  finite  element  methods. 

The  basic  difference  between  the  finite  difference  and  finite  element 
methods  is  that  the  finite  element  method  is  based  upon  approximating  the 
function,  while  the  finite  difference  method  is  founded  upon  approximating  the 
derivatives  of  the  function.  Therefore,  the  finite  difference  method  produces 
solutions  only  at  discrete  points,  while  the  finite  element  method  yields  spatially 
continuous  solutions. 

Also,  the  finite  element  method  offers  numerous  advantages  over  the  finite 
difference  method: 

a.  Anisotropy  and  heterogeneity  of  aquifers  are  easily  taken  care  of. 

b.  Formulation  of  special  formulae  to  incorporate  irregular  boundaries  is 
unnecessary. 

c.  Computer  storage  and  computational  time  can  sometimes  be  saved 
because  often  fewer  nodal  points  are  needed  to  portray  the  region  of 
interest  to  the  same  level  of  accuracy. 
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d.  Irregular  grids  for  handling  different  levels  of  spatial  discretization  in 
different  sections  of  the  region  of  interest  can  be  incorporated. 

e.  The  integral  formulation  used  in  this  method  permits  the  flux  types  of 
the  boundary  conditions  to  come  about  naturally  (Yeh  1987). 

Thus,  the  finite  element  method  is  used  in  this  model.  The  theoretical 
background  as  well  as  numerical  procedures  of  this  method  can  be  found  in  any 
good  finite  element  method  book,  such  as  Istok  (1989),  and  therefore  will  not  be 
described  here.  A  brief  summary  of  the  numerical  procedure  for  applying  the 
finite  element  method  can  be  found  in  Yeh  (1987). 

The  flow  module  of  FEMWATER  includes  three  options  for  solving  the 
finite  element  equations.  In  other  words,  there  are  three  iteration  methods  for 
solving  the  linearized  matrix  equations:  successive  point  iteration,  polynomial 
preconditioned  conjugate  gradient,  and  incomplete  Cholesky  preconditioned 
conjugate  gradient.  Direct  elimination  methods  are  not  used  in  this  report 
because  they  are  impractical  in  dealing  with  large  three-dimensional  problems. 
Because  the  Newton-Raphson  method  will  yield  a  nonsymmetric  matrix,  the 
Picard  method  is  used  to  linearize  the  matrix  equation. 

To  handle  a  large  variety  of  possible  problems,  the  flow  module  for 
FEMWATER  contains  12  optional  numerical  schemes.  Specifically,  the  mixture 
of  schemes  includes  the  combinations  of  (a)  the  three  options  for  solving  the 
resulting  matrix  equation  as  mentioned  in  the  previous  paragraph,  (b)  two 
options  (lumping  and  consistent)  for  handling  the  mass  matrix  resulting  from  the 
storage  term,  and  (c)  two  options  (time-weighted  difference  and  middifference) 
for  approximating  the  time  derivatives.  The  theoretical  background  for  (b)  and 
(c)  may  also  be  found  in  any  respectable  matrix  computation  book  and  in  Yeh 
(1991). 

The  transport  module  for  FEMWATER  also  includes  these  16  options,  plus 
more.  While  the  Galerkin  finite  element  method  is  used  in  the  flow  module,  an 
option  of  two  conventional  finite  element  methods  (Galerkin  and  upstream 
weighting)  or  the  alternative  option  of  a  hybrid  Lagrangian-Eulerian  finite 
element  method  is  provided  in  the  transport  module.  The  main  difference 
between  the  two  conventional  finite  element  methods  is  that  while  the  Galerkin 
finite  element  method  uses  the  base  functions  as  the  weighting  functions,  the 
upstream  weighting  finite  element  method  uses  weighting  functions  different 
from  the  base  functions.  The  advantages  of  using  the  upstream  weighting  finite 
element  method  over  the  Galerkin  finite  element  method  become  apparent  when 
the  advection  terms  in  the  transport  governing  equation  are  equally  important  to 
the  dispersion  terms  (Yeh  and  Ward  1981).  More  details  of  the  two 
conventional  finite  element  methods  may  be  found  in  Yeh  and  Ward  (1981). 


Numerical  Approximation  of  the  Flow  Equations 


98 


Appendix  B  Numerical  Formulation 


Spatial  discretization  with  the  Galerkin  finite  element  method 

When  using  the  finite  element  method,  the  referenced  pressure  head  is 
approximated  by: 

N 

h«h  =  5^hj(t)Nj(x,y,z)  (Bl) 

j=i 


where 


hj  =  amplitude  of  h  at  nodal  point  j 
Nj  =  base  function  at  nodal  point  j 
N  =  total  number  of  nodes 

After  defining  a  residual  and  forcing  the  weighted  residual  to  zero,  the  flow 
equation.  Equation  A28,  is  approximated  as: 


[Nj— FNjdR 

i  Po 


dh; 


dt 


J(VNi)-K-(VNj)dR 


=  J  Ni  ^qdR  -J  (VNi )  •  K  -  -2-  VzdR 

R  Po  R  Po 


(B2) 


Jn-K 


P  1 

Vh  +— Vz  NjdB 
I  Po  j 


In  matrix  form.  Equation  B2  is  written  as: 

[M]|^|  +  [S]{h}  =  {q}  +  {G}  +  {B}  (B3) 

where 

{dh/dt}  =  column  vectors  containing  the  values  of  dh/dt  at  all  nodes 

{/i}  =  column  vectors  containing  the  values  of  h  at  all  nodes 

[M\  =  mass  matrix  resulting  from  the  storage  term 

[5]  =  stiff  matrix  resulting  from  the  action  of  conductivity 

{Q}  =  load  vectors  from  the  internal  source/sink 
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{G}  =  load  vectors  from  the  gravity  force 
{5}  =  load  vectors  from  the  boundary  conditions 
Furthermore,  the  mass  matrix,  [M],  and  stiff  matrix,  [5],  are  described  as: 


Mi,=  SjN;;fEN;dR 

r  0 


eeM,  R 


(B4) 


and 


S«=I  J(VN')  K.(VN;)dR  (B5) 

eeM,  R^ 

where  Re  is  the  region  of  element  e.  Me  is  the  set  of  elements  that  have  a  local 
side  a-P  coinciding  with  the  global  side  i-j,  and  Na  is  the  a-th  local  base 
function  of  element  e. 

In  addition,  the  three  load  vectors,  {Q},  {G},  and  {5},  are  described  as: 

Q.  =  2;jN*^qdR  (B6) 

eeMe  R  Po 


G,=-S  j(VN;).K.;fvzdR  (B7) 

esMt  R,  Po 


and 


-I  jN'n 

B, 


(B8) 


where  Nse  is  the  set  of  boundary  segments  that  have  a  local  node  a  coinciding 
with  the  global  node  i,  and  Be  is  the  length  of  boundary  segment  e. 

In  most  finite  element  work,  the  Darcy  velocity  components  given  in 
Equation  A32  are  calculated  numerically  by  taking  the  derivatives  of  the 
simulated  h  as 


V  = 


(B9) 
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This  formulation  results  in  a  velocity  field  that  is  not  continuous  at  element 
boundaries  and  nodal  points  if  the  variation  of  h  is  other  than  linear  or  constant. 
The  alternative  approach  would  be  to  apply  the  Galerkin  finite  element  method 
to  Equation  A32;  thus  one  obtains 


[t]{vJ={dJ 

(BIO) 

mK}={d,} 

(BID 

[t]{vJ={dJ 

(B12) 

where  the  matrix  [7]  and  the  load  vectors  {Dx],  {Dy},  and  {D^}  are  given  by 


T„=SjNlN;dR 

esMj  R_. 

D„=-S  jN'i-K.Kvh  +  VzjdB 

esM.  IP  J 

jNU  K-j^Vh  +  VzjdB 

eeM.  R^  L  P  J 

Da  =-E  jNJk.K-I^Vh  +  VzidB 

eeM,  R^  IP  J 


(B13) 


(B14) 


(B15) 


(B16) 


where  I4,  Vy,  and  are  the  Darcy  velocity  components  along  the  x-,  y-,  and  z- 
directions,  respectively;  and  i,j,  and  k  are  the  unit  vectors  along  the  x-,  y-,  and  z- 
coordinates,  respectively. 

The  reduction  of  the  partial  differential  equation  (Equation  A28)  to  the  set  of 
ordinary  differential  equations  (Equation  B3)  simplifies  the  evaluation  of 
integrals  on  the  right  side  of  Equations  B4  through  B8  for  every  element  for 
boundary  surface  e.  The  major  tasks  that  remain  to  be  done  are  the  specification 
of  base  and  weighting  functions  and  the  performance  of  integration  to  yield  the 
element  matrices.  Linear  hexahedral  elements  are  employed  in  this 
documentation. 
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Base  and  weighting  functions 

The  construction  of  base  functions  for  hexahedral  elements  is  best 
accomplished  using  the  local  coordinates  (x,  h,  z).  In  the  local  coordinates,  the 
origin^  hexahedral  element  is  mapped  into  a  cubic  whose  comers  are  located  at 
x  =  ±\,h  =  ±1,  and  z  =  ±1  as  shown  in  Figure  Bl. 

For  trilinear  hexahedral  elements,  the  eight  base  functions  are  obtained  by 
taking  the  tensor  product  of  the  three  base  functions  of  the  linear  line  elements  as 

NifeTl,C)  =  -^(l  +  ^^i)(l  +  TlTli)(l  +  CCi)>  i  =  l,2,...,8  (B17) 

Because  the  Galerkin  finite  element  method  is  used  to  solve  the  flow  equations, 
the  set  of  eight  weighting  functions  is  taken  as  the  same  set  of  eight  base 
functions. 


Numerical  integration 

To  complete  the  reduction  of  the  partial  differential  equation  (Equation  B3) 
to  the  ordinary  differential  equation  (Equation  B3),  one  has  to  evaluate  the 
integrals  on  the  right  sides  of  Equations  B4-B8  for  every  element  to  yield  the 
element  mass  matrix  [Af  ]  and  the  stiff  element  matrix  [5*]  as  well  as  the  element 
gravity  colunm  vector  {G*},  the  source/sink  column  vector  {0*}.  and  the 
boundary  column  vector  {B*}  as 


M:p  =  jN^^FN;dR  (B18) 

R,  Po 

Si  =  |(VN;)K(VN;)dR  (B19) 

Re 

* 

Q^  =  jN^^qdR  (B20) 

R,  Po 

G'  =-|(VN;)K.-^VzdR  (B2I) 

R,  Po 


and 
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Since  Equations  B 18-22  are  written  in  global  coordinates  and  the  base 
functions  are  defined  in  local  coordinates,  a  transformation  between  the  global 
and  local  coordinates  is  needed.  The  required  transformation  is  obtained  via  the 
base  functions  as 

8 

x  =  21xjNj(^,Ti,C)  (B23) 

j=i 


yjNjfe-n,;) 


(B24) 


z 


8 


=  SzjNj(tTl,C), 


(B25) 


Because  the  coordinate  transformation  uses  the  base  functions,  the  element  is 
termed  an  “isoparametric”  element. 
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Using  the  transformation  in  Equations  B23-25,  differentiation  of  the  base 
function  with  respect  to  the  global  coordinate  can  be  changed  to  differentiation 
with  respect  to  the  local  coordinate  by 


'9n/ 

aNi' 

dx 

dN, 

■=[jr- 

aN. 

dy 

ari 

dN, 

aN; 

dz 

V.  7 

<  ^  J 

ax 

ay 

az 

ax 

az 

ari 

aij 

ari 

ax 

ay 

az 

ac 

[j]  ’  =  Inverse  of  [J] 


(B26) 


where  [7]  is  the  Jacobian  of  the  transformation.  In  the  mean  time,  a  differential 
volume  written  in  the  local  coordinate  becomes 


1  1  1 

\dR  =  \\\jd^dndi; 

e  -1-1-1 


(B27) 


With  Equations  B26  and  B27,  all  the  integrals  in  Equations  B 18-21  can  be 
reduced  to  the  following  form 

1  1  1 

J  J  J  f(^^nX)jd^dr]d^  (B28) 

-1-1-1 


the  integration  of  which  can  easily  be  carried  out  with  a2x2x2  =  8  point 
Gaussian  quadrature.  The  surface  integration  of  Equation  B22  in  three- 
dimensional  space  is  not  as  straightforward  as  in  two-dimensional  space.  This 
integration  requires  further  elaboration.  Any  surface  integral  of  a  continuous 
function  F(x,y,z)  specified  on  the  surface  S  (Figure  B2)  can  be  reduced  to  the 
area  integration.  Let  I  represent  the  surface  integral: 

I  =  jF(x,y,z)dS  (B29) 

s 


where  the  surface  S  is  given  by  the  following  equation: 


104 


Appendix  B  Numerical  Formulation 


Figure  B2.  A  surface  area  and  its  imbedded  local  coordinates 

z  =  f(x,y)  (B30) 

Let  P  be  any  point  on  the  surface  S  with  coordinates  (x,y,z)  or  (x,  h)  (Figure  B2). 
Then  the  vector  r  from  O  to  P  is  given  by 

r  =  xi  +  yj  +  zk  (B31) 

The  tangent  vectors  to  the  coordinate  curves  x  =  Xo  and  h  =  hg  on  the  surface  S 
are  Sr/Si  and  dr/Sx,  respectively.  The  area  dS  is  given  by: 

9r  9r  „ 

dS=— X— d^dTi  (B32) 

oq  oTi 

where  x  represents  vector  multiplication.  But 
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i  j  k 

dx  dy  dz 

~dl  Jl  Jl 

dx  By  Bz 
Br]  Br\  Brj 

so  that 

ds=^j;+j;+j’d^dii 

where 


3x 

dy 

9x 

dz 

dz 

dy 

Jz  = 

9x 

9x 

dz 

3z 

3ti 

9t| 

dr\ 

9ti 

9ti 

9ti 

Substituting  Equation  B34  into  Equation  B29  gives 
1  1 

j  F(x,y,z)dS  =  J  J  (t)(^,Ti),Jjf+j[+j[d^dTi 

S  -1-1 


Br  Br 


(B33) 


(B34) 


(B35) 


(B36) 


where 

"n)  =  'n).y(t  'n).z(^,  ti))  (B37) 

The  surface  integral  in  Equation  B36  can  easily  be  computed  by  Gaussian 
quadrature. 


Mass  lumping  option 

Referring  to  [M\,  one  may  recall  that  this  is  a  unit  matrix  if  the  finite 
difference  formulation  is  used  in  spatial  discretization.  Hence,  by  proper 
scaling,  the  mass  matrix  can  be  reduced  to  the  finite-difference  equivalent  by 
lumping  (Clough  1971).  Mass  lumping  typically  gives  better  stability  but  less 
accuracy  than  no  lumping.  However,  lumping  can  give  more  accurate  and  stable 
solutions  if  it  is  used  in  conjunction  with  the  central  or  backward-difference  time 
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marching  (Yeh  and  Ward  1980).  Therefore,  options  are  provided  for  the 
lumping  of  the  matrix  [M].  More  explicitly,  [M\  will  be  lumped  according  to: 


M,=  2|SjN:^FN 

r  0 


e6M,\^B=l  R, 


;dR 


if  j  =  1, 


(B38) 


and 


Mij=0  ifj^^l 


(B39) 


Finite  difference  approximation  in  time 

Next,  a  matrix  equation  is  derived  by  integrating  Equation  B3.  For  the  time 
integration  of  Equation  B3,  the  load  vector  {B}  will  be  ignored.  This  load  vector 
will  be  discussed  in  the  next  section  on  the  numerical  implementation  of 
boundary  conditions.  An  important  advantage  of  finite  element  approximations 
over  the  finite  difference  approximations  is  the  inherent  ability  to  handle 
complex  boundaries  and  obtain  the  normal  derivatives  therein.  In  the  time 
dimension,  such  advantages  are  not  evident.  Thus,  finite  difference  methods  are 
typically  used  in  the  approximation  of  the  time  derivative.  Two  time-marching 
methods  are  adopted  in  the  present  flow  model. 

The  first  one  is  the  time-weighted  method  written  as: 

-{h},)+(o[S]{h},.„  +{l-(o)[S]{h}. 

=  {q}+{G} 


where  [A/],  [5],  {Q},  and  {G}  are  evaluated  at  (t  +  wAt).  In  the  Crank-Nicolson 
centered-in-time  approach,  w  =  0.5;  in  the  backward-difference  (implicit 
difference),  w  =  1.0;  and  in  the  forward-difference  (explicit  scheme),  w  =  0.0. 
The  central-Nicolson  algorithm  has  a  trancation  error  of  O(Ai^),  but  its 
propagation-of-error  characteristics  frequently  lead  to  oscillatory  nonlinear 
instability.  Both  the  backward-difference  and  forward-difference  have  a 
truncation  error  of  0(At).  The  backward-difference  is  quite  resistant  to 
oscillatory  nonlinear  instability.  On  the  other  hand,  the  forward  difference  is 
only  conditionally  stable  even  for  linear  problems,  not  to  mention  nonlinear 
problems. 

In  the  second  method,  the  values  of  unknown  variables  are  assumed  to  vary 
linearly  with  time  during  the  time  interval.  At.  In  this  middifference  method,  the 
recurrence  formula  is  written  as: 
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(B41) 


f 

V 


^[M]  +  [S)]{h} 
At  J 


t+At/2 


-^[M]{h},={Q}  +  {G} 


and 


t+At  ~  2{h}  {h} ,  , 


(B42) 


where  [M],  [5],  and  {Q}  are  evaluated  at  (t  +  At/2). 

Equations  B40  and  B41  can  be  written  as  a  matrix  equation 


[T]{h}  =  {Y}, 


(B43) 


where  [7]  is  the  matrix,  {/z}  is  the  unknown  vector  to  be  found  and  represents  the 
values  of  discretized  pressure  field  at  a  new  time,  and  {7}  is  the  load  vector. 
Take,  for  example.  Equation  B40  with  w  =  1.0.  [7]  and  { 7}  represent  the 
following: 


[T]  = 


[M] 

At 


+  [S]. 


(B44) 


and 

{Y}  =  M{h}^+{Q}  +  {G} 


(B45) 


where  {/i}  is  the  vector  of  the  discretized  pressure  field  at  a  previous  time. 


Numerical  implementation  of  boundary  conditions 

The  following  steps  are  the  incorporation  of  boundary  conditions  into  matrix 
equations  by  the  finite  element  method.  For  the  Cauchy  boundary  condition 
given  by  Equation  A40,  simply  substitute  Equation  A40  into  Equation  B22  to 
yield  a  boundary-element  column  vector  {B/}  for  a  Cauchy  segment: 


(B46) 


where  {9/}  is  the  Cauchy  boundary  flux  vector  given  by 
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(B47) 


Bj  Pw 

The  Cauchy  boundary  flux  vector  represents  the  normal  fluxes  through  the  two 
nodal  points  of  the  segment  Be  on  Be-  For  the  Neumann  boundary  condition 
given  by  Equation  A39,  substitute  Equation  A39  into  Equation  B22  to  yield  a 
boundary-element  column  vector  {Bn}  for  a  Neumann  segment; 


(B48) 


where  [qna‘}  is  the  Neumann  boundary  flux  vector  given  by: 


n-K-^^Vz-NX 

Pw  J 


dB; 


a  =  1, 


,4 


(B49) 


which  is  independent  of  pressure  head. 

The  implementation  of  variable-type  boundary  conditions  is  more  involved. 
During  the  iteration  of  boundary  conditions  on  the  variable  boundary,  one  of 
Equations  A41  through  A45  is  used  at  a  node.  If  either  Equation  A42  or  A45  is 
used,  substitute  it  into  Equation  B22  to  yield  a  boundary  element  column  vector 
{B/}  for  a  variable  boundary  segment: 


(B50) 


where  is  the  variable  boundary  flux  given  by: 
C=-jN:.^q,dB,  or 

B,  Pw 

q..=-|N'^q.dB;  a  =  l....,4 

B,  Pw 


(B51) 


Assembling  over  all  Neumann,  Cauchy,  and  variable  boundary  segments  yields 
the  global  boundary  column  vector  {B}  as: 


{B}  =  {q} 


(B52) 


in  which 
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(B53) 


{q}=I{q:}+I{q:}+S{q;} 


eeN„j  eeN„  eeN„ 


where  N„e,  Nee,  and  N^e  are  the  number  of  Neumann  boundary  segments,  Cauchy 
boundary  segments,  and  variable  boundary  segments  with  flux  conditions 
imposed  on  them,  respectively.  The  boundary  flux  {5}  given  by  Equations  B52 
and  B53  should  be  added  to  the  right  side  of  Equation  B43. 

For  the  river  boundary  condition  given  by  Equation  A46,  simply  substitute 
Equation  A46  into  Equation  B22  to  yield  a  boundary-element  column  vector 
{B/}  and  boundary  matrix  [i?/]  for  a  river  boundary  segment: 

{BR}  =  {qr}  and  [B*]  =  [b:]  (B54) 


where  and  {fe/}  are  the  contribution  of  the  Cauchy  boundary  to  the  right 
and  left  sides  of  the  matrix  equation,  respectively.  They  are  given  by 


q 


e 

rot 


hjjdB,  a  =  1,2  and 

Kr  , 

-^NJdB 


(B55) 


At  nodes  where  Dirichlet  boundary  conditions  are  applied,  an  identity 
equation  is  generated  for  each  node  and  included  in  the  matrices  of  Equation 
B43.  The  Dirichlet  nodes  include  the  nodes  on  the  Dirichlet  boundary  and  the 
nodes  on  the  variable  boundary  to  which  either  Equation  A41,  A43,  or  A44  is 
applied. 

After  time  discretization  of  Equation  B3  and  incorporation  of  boundary 
conditions,  the  following  matrix  equation  is  obtained 

[C]{/t}  =  {/?}  (B56) 


where  [C]  is  the  coefficient  matrix  and  {/?}  is  the  known  vector  of  the  right  side. 
For  the  saturated-unsaturated  flow  simulation,  [C]  is  a  highly  nonlinear  function 
of  the  pressure  head  {h}. 


Solution  of  the  matrix  equations 

Equation  B56  is  in  general  a  banded  sparse  matrix  equation.  It  may  be 
solved  numerically  by  either  direct  or  iteration  methods.  In  direct  methods,  a 
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single  solution  operation  sequence  is  performed.  This  would  result  in  an  exact 
solution  except  for  round-off  error.  In  this  method,  one  is  concerned  with  the 
efficiency  and  magnitude  of  round-off  error  associated  with  the  sequence  of 
operations.  On  the  other  hand,  in  an  iterative  method,  one  attempts  to  find  a 
solution  by  a  process  of  successive  approximations.  This  involves  making  an 
initial  guess,  then  improving  the  guess  by  some  iterative  process  until  an  error 
criterion  is  obtained.  Therefore,  in  this  technique,  one  must  be  concerned  with 
convergence,  and  the  rate  of  convergence.  The  round-off  errors  tend  to  be  self- 
corrected. 

For  practical  purposes,  the  advantages  of  direct  methods  are  as  follows:  (a) 
the  efficient  computation  when  the  bandwidth  of  the  matrix  [C]  is  small,  and  (b) 
the  fact  that  no  convergence  problem  is  encountered  when  the  matrix  equation  is 
linear  or  small  convergence  problems  when  the  mass  equation  is  nonlinear.  The 
disadvantages  of  direct  methods  are  the  excessive  requirements  on  CPU  storage 
and  CPU  time  when  a  large  number  of  nodes  is  needed  for  discretization.  On  the 
other  hand,  the  advantages  of  iterative  methods  are  the  efficiencies  in  terms  of 
CPU  storage  and  CPU  time  when  large  problems  are  encountered.  Their 
disadvantages  are  the  requirements  that  the  matrix  [C]  must  be  well  conditioned 
to  guarantee  a  convergent  solution.  For  three-dimensional  problems,  the 
bandwidth  of  the  matrix  is  usually  large;  thus  the  direction  solution  method  is 
not  practical.  Only  iterative  methods  are  implemented  in  FEMWATER.  Four 
iteration  methods  are  used  in  solving  the  linearized  matrix  equation:  (a)  block 
iteration,  (b)  successive  point  iteration,  (c)  polynomial  preconditioned  conjugate 
gradient  method,  and  (d)  incomplete  Cholesky  preconditioned  conjugate  gradient 
method. 

The  matrix  equation.  Equation  B56,  is  nonlinear  because  both  the  hydraulic 
conductivity  and  the  water  capacity  are  functions  of  the  pressure  head  h.  To 
solve  the  nonlinear  matrix  equation,  two  approaches  can  be  taken:  (a)  the  Picard 
method  and  (b)  the  Newton-Raphson  method.  The  Newton-Raphson  method  has 
a  second  order  of  convergent  rate  and  is  very  robust.  However,  the  Newton- 
Raphson  method  would  destroy  the  symmetrical  property  of  the  coefficient 
matrix  resulting  from  the  finite  element  approximation.  As  a  result  the  solution 
of  the  linearized  matrix  equation  requires  extra  care.  Many  of  the  iterative 
methods  will  not  give  a  convergent  solution  for  the  nonsynunetric  linearized 
matrix  equation.  Thus,  the  Picard  method  is  used  in  this  report  to  solve  the 
nonlinear  problems. 

In  the  Picard  method,  an  initial  estimate  is  made  of  the  unknown  {/z}.  Using 
this  estimate,  the  coefficient  matrix  [C]  is  then  computed  and  the  linearized 
matrix  equation  solved  using  linear  algebra.  The  new  estimate  is  now  obtained 
by  the  weighted  average  of  the  new  solution  and  the  previous  estimate: 

=  (B57) 


where 
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=  new  estimate 
{/i*}  =  previous  estimate 
{/i}  =  new  solution 
(0  =  iteration  parameter 

The  procedure  is  repeated  until  the  new  solution  {/z}  is  within  a  tolerance  error. 
If  w  is  greater  than  or  equal  to  0  but  is  less  than  1,  the  iteration  is 
underrelaxation,  ff  w  =  1,  the  method  is  the  exact  relaxation.  If  w  is  greater  than 
1  but  less  than  or  equal  to  2,  the  iteration  is  termed  overrelaxation.  The 
underrelaxation  should  be  used  to  overcome  cases  of  non-convergence  or  slow 
convergence  rates  due  to  fluctuations  rather  than  due  to  “divergent” 
computations.  Overrelaxation  should  be  used  to  speed  up  the  convergence  rate 
when  it  decreases  monotonically. 

In  summary,  there  are  sixteen  optional  numerical  schemes  here  to  deal  with  a 
wide  range  of  problems.  These  are  the  combinations  of  (a)  two  ways  of  treating 
the  mass  matrix  (lumping  and  no-lumping);  (b)  two  ways  of  approximating  the 
time  derivatives  (time-weighting  and  middifference),  and  (c)  four  ways  of 
solving  the  linearized  matrix  equation. 


Transport  Equation 


Spatial  discretization  with  the  weighted  residual  finite  element 
method 

Let  Cj  be  approximated  by  a  finite  element  interpolation  as 

N 

C  =  C  =  2  C  j  (t)N  j  (x,  y,  z)  (B58) 

j=i 


Substituting  Equation  B58  into  Equations  A57,  A61,  and  A63,  and  forcing  a 
weighted  residual  to  zero  gives  the  following  ordinary  differential  equations: 

a.  For  the  conventional  finite  element  approach: 
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|^}+  lw,\-VN,dR  {c}- 


JVA?,  flD  VWjdS  {C}+|jw,  ^eK,+p,K,^'j  JV,<«|{C}' 

{p,[(a'f+A)(e+ftf)]w,.«}{C}+ 

■  +^'e)§-^v.vW|{c}= 

R  [P  ^  P  VPoJj] 

Jw,[^-a  +  i:,)p.  5-^Cj<iR  + 

J  N.qCJR  +  J  iV,.  n  •  00  •  VCdB 

R  B 

b.  For  the  Lagrangian-Eulerian  approach  with  a  linear  isotherm: 


]Nj^e  +  p,^N,dR  1^1+  JVAf,.®  VJV;<iR  {C}+ 
jjiV,  ieK^+p,K,^'j  Ar,<fl!|{C}  + 

[i  I  P  ^  P  VPoJjj 

.  r  ^  r  f 

[A^,.  -a  +  KJp,  S——C  dR  +  \N,qCJR+\N,n-eD-VCdB 

i  I  y  dc  jj  i 

c.  For  the  Lagrangian-Eulerian  approach  with  a  nonlinear  isotherm: 


(B59) 


(B60) 
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J  N,qC.JR  +  J  A^,.  n  eD  •  ^CdB 

R  B 


Equations  B59-B61  are  written  in  matrix  form  as: 

a.  For  the  conventional  finite  element  approach: 

[M]{^}  +  Ca] + [D]  +  [K]XC}  =  {q}  +  {B}  (B62) 

b.  For  the  Lagrangian-Eulerian  approach  with  the  linear  isotherm: 

[M]|^^}  +  ([D]  +  [K]){C}  =  {q}  +  {B}  (B63) 

c.  For  the  Lagrangian-Eulerian  approach  with  nonlinear  isotherms: 

[m,  ]|^|  +  [m,  ]|^|  j  +  ([D] + [KJXc}  =  {Q} + {B}  (B64) 

where 

{C}  =  vector  whose  components  are  the  concentrations  at  all  nodes 

{dC/dt}  =  derivative  of  {C}  with  respect  to  time 

[M],  [M;]  =  mass  matrices  associated  with  the  material  derivative  term 
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[M2]  =  mass  matrix  associated  with  the  partial  derivative  term 
[D]  =  stiff  matrix  associated  with  the  dispersion  term 
[A]  =  stiff  matrix  associated  with  the  advection  term 
[/n  =  stiff  matrix  associated  with  all  the  &st-order  terms 
{Q}  =  load  vector  associated  with  all  zero-order  derivative  terms 
{5}  =  load  vector  associated  with  boundary  conditions 
These  matrices  and  vectors  are  given  as: 


M,j=S|N;fe  +  Pi^KdR  (B65) 

eeM  V  / 

Muj=SjN^eN;dR  (B66) 

eeMR^ 

M„=SjN'fpi^W  (B67) 

Aij  =EjN:V-VN;dR  (B68) 

eeMR^ 

Dij  =  S  J •  VN;dR  (B69) 

eeMR^ 
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eeM  R 


dh  I 

\k 

-{X  +  a—  +  K^)p,\ 

1  dC  J 

dR  (B71) 


B;  =  -X  J  VN>  •  (-  eD  ■  VC)dB 


(B72) 


csMb 


where  Cm  is  the  concentration  of  the  source. 


Base  and  weighting  functions 

For  the  flow  case,  the  weighting  functions  are  taken  as  the  same  set  as  the 
base  functions.  However,  in  the  transport  formulation  using  an  Eulerian  finite 
element  approach,  sometimes  it  is  advantageous  to  use  weighting  functions 
which  are  one  or  two  orders  higher  than  the  base  functions:  (N+1)  or  (N+2) 
upstream  weighting.  This  section  will  present  only  the  N+1  upstream  weighting 
functions.  Recently,  the  N+2  weighting  functions  have  been  the  subject  of 
several  investigations.  The  success  of  the  N+2  weighting  is  still  under 
investigation.  They  will  not  be  included  here.  First  define,  for  the  line  element, 
the  following  N+1  upstream  weighting  functions 


F,(4,a)  =  N.(|)-a|(l+#-4) 

(B73) 

F,(ia)  =  N,{4)  +  a|(l  +  4)(l-5) 

(B74) 

where  a  is  the  weighting  factor  along  the  line  from  node  1  to  node  2  (Figure  B3). 


Figure  B3.  Weighting  factor  along  a  line  element 

Then  the  weighting  functions  are  obtained  by  an  appropriate  tensor  product: 
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W,=F.{ta.)F,(ll,P,)F,{CY,) 

(B75) 

W,=F,(4,a,)F,(ll,P,)F,(C,7,) 

(B76) 

W,=F,(ta,)F,(ll,Pi)F,fcY,) 

(B77) 

W.=F.(ta,)F,(ll,p,)F,(?,7,) 

(B78) 

W,=F,(ta3)F,(ll,P3)F,(?,7,) 

(B79) 

W3=F3(ta3)F,(T|,P3)F3(?,73) 

(B80) 

W,  =F3(ta,)F3(ll,p,)F3fcY,) 

(B81) 

W,=F,fe,a,)F3(ll.p3)F!fcY3) 

(B82) 

where  a’s,  13’ s,  and  /s  are  the  weighting  factors  along  the  side  given  in  Figure 
B4. 


Numerical  integration 


To  reduce  the  partial  differential  equations.  Equations  A57,  A61,  and  A63  to 
ordinary  differential  equations.  Equations  B62-B64,  one  has  to  evaluate  the 
integrals  on  the  right  sides  of  Equations  B65-B72  for  every  element  to  yield  the 
element  mass  matrices  [Af],  [Mi‘],  and  [M^*]  and  the  stiff  element  matrices  [A*], 
[D%  and  [^]  as  well  as  the  source/sink  column  vector  {Q‘}  and  the  boundary 
colunrn  vector  {B*}  as: 


(B83) 


MU  =  j  NjeNJdR 

Rc 


(B84) 
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Figure  B4.  Upstream  weighting  factors  along  12  sides  of  a  hexahedral  element 


.  f  J  dS\  , 

jN'pdR 

(B85) 

a;i,  =  jN^VVN^dR 

(B86) 

Rt 

Drt  =  ^N^eD■VN5ciR 

(B87) 

Rc 


Klf=\K 

R, 


C  dS^ 

Pb^s^ 


V 


dh 

+  {(X~z:~  +  X). 
at 


dCj 
d  +  Pb 


dC 


P*  /  d  _  dh 
+—q-(a —  +  pd)— 
p  at 


Po 


V-V 


\NldR 


(B88) 


Appendix  B  Numerical  Formulation 


<x.  =  \k 


dh 

-i^  +  CC~^  +  ^s)Pb 


dS  ^ 
dC^J 


+  ^C',.n 


dR 


(B89) 


=-jN>(-0DVC)dB 


(B90) 


Following  the  procedures  presented  in  the  section  “Numerical  integration,” 
Equations  B83-B90  are  first  transformed  in  terms  of  local  coordinates.  Then  the 
resulting  equations  are  integrated  with  Gaussian  quadrature.  The  transformation 
between  the  global  and  local  coordinates  is  given  by  Equations  B23-B25 
resulting  in  isoparametric  elements.  The  surface  integration  from  the  boundary 
conditions  also  follows  that  presented  in  “Numerical  integration.” 


Mass  lumping  option 

As  with  the  solution  of  the  flow  equations,  a  consistent  mass  matrix  or  mass 
lumping  option  can  be  used  when  the  Eulerian  approach  is  used.  Although  a 
consistent  mass  matrix  option  can  be  used  when  the  hybrid  Lagrangian-Eulerian 
approach  is  taken,  a  mass  lumping  scheme  is  more  appropriate  and  easier  to 
implement. 


Finite  difference  approximation  in  time 

When  the  Eulerian  approach  is  taken  in  approximating  the  governing 
equations,  either  time-weighted  difference  or  middifference  can  be  used  as  in  the 
section  “Finite  difference  approximation  in  time.”  However,  when  the 
Lagrangian-Eulerian  approach  is  taken,  the  time  integration  is  different  from  the 
flow  problem.  Although  there  is  still  a  choice  of  time-weighted  difference  or 
middifference,  the  time-weighted  difference  scheme  is  preferred.  In  the 
following,  the  time  integration  for  the  Lagrangian-Eulerian  approach  is 
demonstrated.  As  in  the  time  integration  of  the  flow  equations,  the  boundary 
load  vector  will  be  ignored  in  the  time  integration  of  the  transport  equations  in 
this  section.  This  load  vector  will  be  discussed  in  the  next  section. 

In  the  Lagrangian-Eulerian  approach.  Equations  B63  and  B64  are  integrated 
along  the  characteristic  lines.  First  Equation  B63  will  be  integrated  for  the  linear 
isotherm  case.  Then  Equation  B64  will  be  integrated  for  nonlinear  isotherms. 
The  time-weighted  integration  of  Equation  B63  yields 

M({c-  }-  {c-})+ (o([D]  +  [K]){C"*‘ }  + 

{i-co)([d]+[k]){c'}={q} 
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where  -4Tis  the  time-step  size  (the  determination  of  ^Twill  be  explained  later), 

{ is  the  concentration  vector  containing  the  concentration  at  all  nodes  at  the 
new  time  n+1,  and  {C*}  is  the  Lagrangian  concentration  vector.  The  Lagrangian 
concentration  {C*}  is  computed  by  the  backward  method  of  characteristics  as 
follows. 


t+At 


X*  =Xi-  Jv,dt 

t 

c-=i;cj(t)Nj(x:) 


(B92) 


where  x,  *  (the  Lagrangian  point)  is  the  location  of  a  fictitious  particle  originating 
at  time  t  and  arriving  at  the  node  x,  at  time  t  +  At.  Ci(t)  is  the  concentration  value 
at  node  j  at  time  t  and  N/x,*)  is  the  interpolation  function  associated  with  nodej 
evaluated  at  the  Lagrangian  point  x,  *.  If  x,  *  is  located  in  the  interior  of  the 
region  of  interest.  Atm  Equation  B91  is  defined  as 


At  =  At 


(B93) 


If  Ax*  is  located  outside  the  region  of  interest,  a  Atfx,  *)  as  a  function  of  x,  *  must 
be  found  such  that 


Xi  =Xi 


t+AT(xi) 

Jv,dt 


(B94) 


will  be  located  on  the  boundary.  Thus,  it  can  be  seen  that  At  is  less  than  or  equal 
to  At. 


For  the  nonlinear  isotherm  case.  Equation  B64  can  be  integrated  to  yield 


[m.]  ^  [mJ 


At  At 


{C“^'}  +  0)([D]  +  [K]){C°^'} 


1m.],  ,  [mJ 

At  ^  ^  At 


{Cn} 


(B95) 


(l-a))([D]-h[K]){c*}  +  {Q} 


The  computation  of  Ar  and  the  Lagrangian  concentrations  C*  in  Equation  B95 
follows  Equations  B92-B94  but  with  replaced  by  Vf. 
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Numerical  implementation  of  boundary  conditions 


To  incorporate  the  boundary  conditions,  the  right  side  of  Equation  B90  must 
be  evaluated  for  every  boundary  segment  Be  to  yield  the  load  vector  {ff}\ 

=-jN>-(-0D-VC)dB,  a  =  l,...,4  (B96) 

Be 

For  the  Neumann  boundary  condition  given  by  Equation  A67,  simply  substitute 
Equation  A67  into  Equation  B90  to  yield  a  boundary-element  column  vector 
{Bn‘}  for  a  Neumaim  segment; 


(B97) 


where  {9/}  is  the  Neumarm  boundary  flux  vector  given  by 

q»=-|N'qNdB,  a  =  l,...,4  (B98) 

Be 

This  Neumann  boundary  flux  vector  represents  the  normal  fluxes  through  the 
two  nodal  points  of  the  segment  Be  on  B„. 

For  the  Cauchy  boundary  condition  given  by  Equation  A68,  Equation  B90 
may  be  rewritten  in  the  following  form: 

=-jN>-(VC-eD-VC)dB-i-jN>-VCdB,  a  =  l,...,4  (B99) 

Be  B, 


The  concentration  on  the  boundary  segment  Be  can  be  approximated  by 


(BlOO) 


Substituting  Equations  A68  and  BlOO  into  Equation  B99  gives  a  boundary- 
element  column  vector  {B/}  for  a  Cauchy  segment: 

{B:}={q:}+[v']{C}  (BlOl) 

in  which  the  Cauchy  boundary  flux  vector  {5/}  and  the  Cauchy  boundary  matrix 
[Vf]  from  the  normal  velocity  component  are  given  by 
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(B102) 


{q»}  =  -jN'q.dB,  a=l,....4 
and 

V:.p  =  jN>-VN;ciB,  a  =  l,...,4  and  P  =  l,...,4 

b. 


Segments  on  which  the  variable  boundary  conditions  are  imposed  are  the 
flow-through  boundaries  on  which  the  flow  direction  is  not  known  a  priori. 

When  the  flow  is  directed  into  the  region,  Cauchy  boundary  conditions  will  be 
used.  The  boundary  element  column  vector  {By}  for  a  variable  boundary 
segment  can  be  obtained  similarly  to  {B/}: 

{B:}={q;}+[v;]{c}  (BIOS) 


in  which  the  variable  boundary  flux  vector  {^v*}  and  the  variable  boundary 
matrix  [V/]  from  the  normal  velocity  component  are  given  by: 

and  (B104) 

v:^=\Kn-yNldB,  a  =  l,...,4  and  j8  =  l,...,4 

b. 


where  C,„  is  the  total  dissolved  concentration  of  the  incoming  fluid.  When  the 
flow  is  directed  out  from  the  region,  both  {q^}  and  [FvT  are  set  equal  to  0. 

Assembling  over  all  Neumann,  Cauchy,  and  variable  boundary  segments,  the 
global  boundary  column  vector  {B}  is  obtained  as: 

{B}  =  {q}-H[V]{C}  (BIOS) 


in  which 


{q}=  S{q'.}4E{q:}4l{q'.} 

eeN^  eeN„  eeN,, 

[v]=  ZK'l+SK'} 


eeN„ 


eeN„ 


(B106) 


where  Nm,  N^,  and  Nyg  are  the  number  of  Neumann,  Cauchy,  and  variable 
boundary  segments,  respectively. 
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At  nodes  where  Dirichlet  boundary  conditions  are  applied,  an  identity 
equation  is  generated  for  each  node  and  included  in  the  matrices  of  Equation 
B91  for  the  case  of  linear  isotherms  or  Equation  B95  for  the  case  of  nonlinear 
isotherms.  The  detailed  method  of  applying  this  type  of  boundary  condition  can 
be  found  in  Wang  and  Connor  (1975). 

Boundary  conditions  need  to  be  implemented  in  the  computation  of  the 
Lagrangian  concentrations  {C*}.  Neumann  boundary  conditions  are  normally 
applied  to  the  boundary  when  flow  is  directed  out  from  the  region  of  interest. 

On  the  Neumann  boundary,  the  backtracking  would  locate  x,*  in  the  interior  of 
the  domain;  hence  the  Lagrangian  concentration  at  the  ith  Neumann  boundary 
node  is  simply  computed  via  interpolation.  On  the  Dirichlet  boundary  nodes,  the 
Lagrangian  concentration  is  simply  set  to  the  specified  value. 

On  the  variable  boundary,  boundary  conditions  need  not  be  implemented  if 
the  flow  is  directed  out  from  the  region.  If  the  flow  is  directed  into  the  region, 
the  concentration  of  incoming  fluid  is  specified.  An  intermediate  concentration 
C**  is  calculated  according  to 

C"  =  J  Nj  V„Ci,dB  /Nj  V„dB  ,  (B107) 

b; 

where  Q**  is  the  concentration  due  to  the  boundary  source  at  the  boundary  node 
i,  V„  is  the  normal  vertically  integrated  Darcy’s  velocity,  and  C,„  is  the 
concentration  of  the  incoming  fluid. 

Cauchy  boundary  conditions  are  normally  applied  to  the  boundary  where 
flow  is  directed  into  the  region,  where  the  material  flux  of  incoming  fluid  is 
specified.  The  intermediate  concentration  is  thus  calculated  according  to 

C  =  J  N.q^dB  /  J  NyjB  ,  (B 108) 

Be  B, 


where  C,**  is  the  concentration  due  to  Cauchy  fluxes  at  the  boundary  node  i,  V„ 
is  the  normal  Darcy’s  velocity,  and  qc  is  the  Cauchy  flux  of  the  incoming  fluid. 

The  Lagrangian  concentration  is  obtained  by  using  the  value  C,  **  and  C," 

(the  concentration  at  previous  time-step)  as  follows: 

J  N,eN,CrdB  +  jN,p,K,N,C;dB 

C*  =  - } — ; — - ^ - - for  the  linear  isotherm  (B109) 

jN,(e  +  p,K,)dB 


C*  =  C**  for  the  nonlinear  isotherm 


(BllO) 
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Solution  of  the  matrix  equations 


Because  the  Lagrangian-Eulerian  approach  results  in  a  symmetric  positive 
definite  matrix,  the  system  of  algebraic  equations  can  be  solved  by  any  of  the 
four  options:  block  iteration,  successive  point  iteration,  polynomial 
preconditioned  conjugate  gradient,  and  incomplete  Cholesky  preconditioned 
conjugate  gradient  methods.  For  the  Eulerian  approach,  however,  the  block 
iteration  and  successive  point  iteration  methods  are  the  preferred  choice  for 
solving  the  matrix  equation.  When  the  advection  transport  is  dominant,  the  two 
basic  iteration  methods  with  underrelaxation  are  very  effective  in  reducing  the 
number  of  iterations  required  for  a  convergent  solution  (Yeh  1985). 
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Appendix  C  Output  File 
Formats 


Introduction 

The  output  from  FEMWATER  consists  of  a  printed  output  file  and  a  number 
of  solution  files.  The  printed  output  file  is  a  text  file  listing  a  summary  of  the 
input,  iteration  and  convergence  data,  and  solution  summaries.  The  solution 
files  are  used  for  post-processing  or  as  initial  conditions  for  subsequent  runs. 

The  content  and  format  of  the  solution  files  are  described  in  this  appendix. 

There  are  four  solution  files  that  can  be  output  from  FEMWATER:  pressure 
head,  velocity,  moisture  content  (nodal),  and  concentration.  All  of  the  files  can 
be  output  in  either  binary  or  text  format.  The  options  for  selecting  which  files 
are  to  be  output  and  which  format  is  used  are  specified  on  the  OC4  card  in  the 
model  file  (see  section,  “Save  options  (OC4),”  in  main  text). 


Data  Set  Files 


The  pressure  head,  moisture  content  (nodal),  concentration,  and  velocity 
files  are  all  written  using  the  standard  Department  of  Defense  Groundwater 
Modeling  System  (GMS)  data  set  file  format.  The  data  set  file  formats  are  two 
of  the  standard  file  formats  used  by  the  GMS.  Multiple  data  sets,  including  both 
scalar  and  vector,  can  be  included  in  one  file.  However,  the  data  set  files  used 
by  FEMWATER  for  input  and  output  are  assumed  to  contain  one  data  set  per 
file. 


With  the  scalar  data  sets,  one  scalar  value  is  listed  per  node,  in  sequential 
order  based  on  the  node  ID’s.  If  the  data  set  is  transient,  a  complete  set  of  scalar 
values  is  listed  for  each  time-step.  With  vector  data  sets,  the  x-,  y-,  and  z- 
components  of  the  vector  (velocity  in  this  case)  are  listed  on  a  node  by  node 
basis. 
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Data  set  files  can  be  saved  in  either  text  or  binary  format.  The  binary 
format  results  in  smaller  files  and  can  result  in  large  reductions  in  the  disk  space 
required,  particularly  when  doing  transient  analyses  on  large  meshes. 


Text  Data  Set  File  Format 

A  summary  of  the  text  version  of  the  scalar  data  set  file  format  is  shown  in 
Figure  Cl. 


DATASET 

/*  File  type  identifier  */ 

OBJTYPE  type 

r  Type  of  objects  */ 

BEGSCL 

/*  Beginning  of  scalar  data  set  */ 

ND  numdata 

/*  Number  of  data  values  */ 

NC  numcells 

/*  Number  of  cells  or  elements  */ 

NAME  “name” 

/*  Data  set  name  7 

TS  istat  time 

r  Time-step  of  the  following  data  7 

stat, 

stat^ 

val, 

valj 

/*  Repeat  TS  card  for  each  time-step  */ 
ENDDS 

r  End  of  data  set  7 

BEGVEC 

r  Beginning  of  vector  data  set  7 

ND  numdata 

r  Number  of  data  values  7 

NC  numcells 

/*  Number  of  cells  or  elements  7 

NAME  “name” 

/*  Data  set  name  7 

TS  istat  time 

/*  Time-step  of  the  following  data  7 

stat, 

statj 

staU,^. 

V.,  v„  v„ 

V  V  V 
''x2  V  ’'zs 

V  V  V 

*xnumdBtB  ’ynumdsta  *zrxjmdMi 

/*  Repeat  TS  card  for  each  time-step  */ 
ENDDS 

r  End  of  data  set  7 

/*  Repeat  BEGSCL  and  BEGVEC  sequences  for  each  data  set.  */  I 

Figured.  Text  data  set  file  format 

The  first  line  of  the  file  is  a  card  without  any  fields  that  serves  as  a  file  type 
identifier. 


Card  Type 

DATASET 

Description 

File  type  identifier. 

Required 

YES 
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The  next  line  is  an  identifier  that  tells  GMS  which  of  the  objects  the  data  set 
is  associated  with.  For  FEMWATER,  the  object  type  should  always  be  three- 
dimensional  (3-D)  mesh. 


Card  Type 

OBJTYPE 

Description 

Identifies  the  type  of  objects  that  the  data  sets  in  the  file  are  associated  with. 

Reauired 

YES 

Format 

OBJTYPE  type 

Sample 

OBJTYPE  tin 

Field 

Variable 

Value 

Description 

1 

type 

tin 

TINS. 

mesh2d 

2-D  meshes. 

grid2d 

2-D  grids. 

scat2d 

2-D  scatter  points. 

meshSd 

3-D  meshes. 

gridSd 

3-D  grids. 

scatSd 

3-D  scatter  points. 

To  begin  a  data  set,  either  a  BEGSCL  or  BEGVEC  card  is  required, 
depending  on  the  type  of  data  set. 


Card  Type 

BEGSCL 

Description 

Marks  the  beginning  of  a  set  of  cards  defining  a  scalar  data  set. 

Required 

YES 

Card  Type 

BEGVEC 

Description 

Marks  the  beginning  of  a  set  of  cards  defining  a  yector  data  set. 

Required 

YES 

The  pair  of  lines  are  the  ND  and  NC  cards.  These  cards  are  used  to  specify 
the  number  of  nodes  and  elements  in  the  mesh. 


Card  Type 

ND 

Description 

The  number  of  data  values  that  will  be  listed  per  time-step.  This  number  should 
correspond  to  the  number  of  nodes  for  a  3-D  mesh. 

Required 

YES 

Format 

ND  numdata 

Sample 

ND  10098 

Field 

Variable 

Value 

Description 

1 

numdata 

+ 

The  number  of  items.  At  each  time-step,  numdata 
values  are  listed. 

Card  Type 

NC 

Description 

This  number  should  correspond  to  the  number  of  elements  for  a  3-D  mesh. 

Required 

YES 

Format 

NC  numcells 

Sample 

NC  3982 

Field 

Variable 

Value 

Description 

1 

numcells 

+ 

The  number  of  elements. 
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The  next  line  is  for  the  name  of  the  data  set. 


Card  Type 

NAME 

Description 

The  name  of  the  data  set. 

Required 

YES 

Format 

NAME  “name” 

Sample 

NAME  “Total  head” 

Field 

Variable 

Value 

Description 

1 

“name” 

str 

The  name  of  the  data  set,  in  doubie  quotes. 

A  data  set  can  contain  multiple  solutions,  each  solution  representing  a 
complete  set  of  nodal  values  at  a  particular  time-step.  The  TS  card  is  used  to  list 
a  time  value  and  the  corresponding  set  of  scalar  or  vector  values.  If  the  solution 
is  steady  state,  only  one  TS  card  is  used  and  the  time  value  is  set  to  0.0. ,  If  the 
solution  is  transient,  the  TS  card  is  repeated  once  for  each  time-step. 


Card  Type 

TS  (SCALAR) 

Description 

Defines  a  set  of  scalar  values  associated  with  a  time-step.  Shouid  be  repeated  for 
each  time-step. 

Required 

YES 

Format 

TS  istat  time 

statt 

stat2 

statnumcells 

vail 

val2 

valnumdata 

Sample 

TS  012.5 

34.5 

74.3 

58.4 

72.9 

Field 

Variable 

Value 

Description 

1 

istat 

0 

1 

Use  status  flags  from  previous  time-step.  For  first  time- 
step,  this  value  indicates  that  all  cells  are  active. 

Status  flags  will  be  listed. 

2 

time 

+ 

The  time-step  value.  This  number  is  ignored  if  there  is 
only  one  time-step. 

stat 

0 

1 

Inactive. 

Active. 

One  status  flag  should  be  repeated  per  line  for  each  cell 
or  element.  These  flags  are  included  only  when  ISTAT  = 

1. 

val 

± 

The  scalar  values,  one  per  line. 
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Card  Type 

TS  (VECTOR) 

Description 

Defines  a  set  of  vector  values  associated  with  a  time-step.  Should  be  repeated  for 

each  time-step. 

_ 

Repuired 

YES  I 

Format 

TS  istat  time 

statt 

stat2 

statnumcells 
vxl  vyl  vz1 
vx2  vy2  vz2 

vxn  wn  vzn 

Sample 

TS012.5 

34.5  74.4  634.4 

74.3  643.4  636.3 

58.4  745.4  346.3 

72.9  734.3  345.3 

Field 

Variable 

Value 

Description 

1 

istat 

0 

Use  status  flags  from  previous  time-step.  For  first  time- 
step,  this  value  indicates  that  all  cells  are  active. 

1 

Status  flaas  will  be  listed. 

2 

time 

+ 

The  time-step  value.  This  number  is  ignored  if  there  is 
onlv  one  time-step. 

stat 

0 

Inactive. 

1 

Active. 

One  status  flag  should  be  repeated  per  line  for  each  cell 
or  element.  These  flags  are  included  only  when  ISTAT  = 
1. 

vx  w  vz 

± 

The  vector  values,  one  set  per  line. 

Each  data  set  should  be  terminated  with  an  ENDDS  card. 


Card  Type 

ENDDS 

Description 

Signals  the  end  of  a  set  of  cards  defining  a  data  set. 

Required 

YES 

Binary  Data  Set  File  Format 

Data  set  files  can  also  be  written  in  binary  format.  The  binary  format  is 
patterned  after  the  text  format  in  that  information  is  written  to  the  file  in  groups 
or  “cards.”  There  are  a  few  additional  cards  in  the  binary  version,  but  most  of 
the  cards  are  identical  to  cards  in  the  text  version  except  that  the  card  identifiers 
are  not  written  with  character  strings.  Rather,  they  are  written  with  integer 
codes. 

The  cards  in  the  binary  data  set  file  are  as  follows: 
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Card  Type 

VERSION 

Card  ID 

3000 

Description 

File  type  identifier. 

Required 

Yes 

Card  Type 

OBJTYPE 

Card  ID 

too 

Description 

Identifies  the  type  of  objects  that  the  data  sets  in  the  file  are  associated  with. 

Required 

YES. 

Field 

Variable 

Size 

Value 

Description 

1 

id 

4  byte  int 

1 

TINS. 

2 

Boreholes. 

3 

2-D  meshes. 

4 

2-D  grids. 

5 

2-D  scatter  points. 

6 

3-D  meshes. 

7 

3-D  grids. 

8 

3-D  scatter  points. 

Card  Type 

SFLT 

Card  ID 

110 

Description 

Identifies  the  number  of  bytes  that  will  be  used  in  the  remainder  of  the  file  for  each 
floatinq  point  value  (4,  8,  or  16). 

Required 

YES 

Field 

Variable 

Size 

Value 

Description 

1 

sfit 

4  byte  Int 

+ 

4,  8,  or  16 

Card  Type 

SPLG 

Card  ID 

120 

Description 

Identifies  the  number  of  bytes  that  will  be  used  In  the  remainder  of  the  file  for  status 
flags  (1 , 2,  or  4). 

Required 

YES 

Field 

Variable 

Size 

Value 

Description 

1 

sfig 

4  byte  int 

+ 

1, 2,  or  4 

Card  Type 

BEGSCL 

Card  ID 

130 

Description 

Marks  the  beginning  of  a  set  of  cards  defining  a  scalar  data  set. 

Required 

YES 

Card  Type 

BEGVEC 

Card  ID 

140 

Description 

Marks  the  beginning  of  a  set  of  cards  defining  a  vector  data  set. 

Required 

YES 
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Card  Type 

NUMDATA 

Card  ID 

170 

Description 

The  number  of  data  values  that  will  be  listed  per  time-step.  This  number  should 
corresDond  to  the  number  of  nodes. 

Reauired 

YES 

Reid 

Variable 

Size 

Value 

Description 

1 

numdata 

4  byte  int 

+ 

The  number  of  items.  At  each  time-step, 
numdata  are  listed. 

Card  Type 

NUMCELLS 

Card  ID 

180 _ _ _ 

Description 

This  number  should  correspond  to  the  number  of  elements. 

Reauired 

YES 

Reid 

Variable 

Size 

Value 

Description 

1 

numcells 

4  byte  int 

+ 

The  number  of  elements  or  cells. 

Card  Type 

NAME 

Card  ID 

190 

Description 

The  name  of  the  data  set. 

Reauired 

YES 

Reid 

Variable 

Size 

Value 

Description 

1 

name 

80  bytes 

str 

The  name  of  the  data  set.  Use  one  character 
per  byte.  Mark  the  end  of  the  string  with  the  \0 
character. 

Card  Type 

TS  (SCALAR) 

Card  ID 

200 

Description 

Defines  a  set  of  scalar  or  vector  values  associated  with  a  time-step.  Should  be 
repeated  for  each  time-step. 

Reauired 

YES 

Reid 

Variable 

Size 

Value 

Description 

1 

istat 

SFLG  int 

0 

1 

Use  status  flags  from  previous  time-step.  For 
the  first  time-step,  this  value  indicates  that  all 
cells  are  active. 

Status  flaqs  will  be  listed. 

2 

time 

SFLT  real 

+ 

The  time-step  value.  This  number  is  ignored  if 
there  is  only  one  time-step. 

stat 

SFLG  int 

0 

1 

Inactive. 

Active. 

One  status  flag  should  be  repeated  for  each 
cell  or  element.  These  flags  are  included  only 
when  istat  =  1 . 

val 

SFLT  real 

The  scalar  values.  Repeat  numdata  times. 
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Card  Type 

TS  (VECTOR) 

Card  ID 

200 

Description 

Defines  a  set  of  scalar  or  vector  values  associated  with  a  time-step.  Should  be 

repeated  for  each  time-step. 

Required 

YES 

Field 

Variable 

Size 

Value 

Description 

1 

istat 

SFLG  int 

0 

Use  status  flags  from  previous  time-step.  For 
the  first  time-step,  this  value  indicates  that  all 
cells  are  active. 

1 

Status  flags  will  be  listed. 

2 

time 

SFLT  real 

+ 

The  time-step  value.  This  number  is  ignored  if 
there  is  only  one  time-step. 

stat 

SFLG  int 

0 

Inactive. 

1 

Active. 

One  status  flag  should  be  repeated  for  each 
cell  or  element.  These  flags  are  included  only 
when  istat  =  1 . 

vx  w  vz 

SFLT  real 

The  vector  values.  Repeat  numdata  times. 

Card  Type 

ENDDS 

Card  ID 

210 

Description 

Signals  the  end  of  a  set  of  cards  defining  a  data  set. 

Required 

YES 
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